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ABSTRACT 

The  application  of  analysis  lattice  filters  to  the 
problem  of  determining  the  input  to  a  system  from 
observations  of  the  system's  output  (i.e.,  deconvolution )  is 
discussed.  Both  linear  and  nonlinear  systems  are 
considered.  Lattice  filter  modeling  algorithms  (Levinson 
and  Schur)  are  presented. 

The  theory  of  least-squares  inverse  filters  is  reviewed. 
This  leads  to  a  discussion  of  the  lattice  filter,  which  in 
turn  leads  to  the  Generalized  Lattice  Theory.  The 
Generalized  Lattice  Theory  is  then  used  to  develop  a 
nonlinear  lattice  structure.  Simulations  show  that  the 
nonlinear  lattice  is  an  effective  inverse  filter  for  both 
linear  and   nonlinear  systems. 


TABLE  OF  CONTENTS 

I.   INTRODUCTION  7 

II.   LINEAR  LEAST-SQUARES  DECONVOLUTION  TECHNIQUES-  10 

A.  INTRODUCTION 10 

B.  KALMAN  FILTER 14 

C.  LEAST-SQUARES  INVERSE  FILTER  26 

D.  LATTICE  FILTER 38 

1  .   Introduction 38 

2.  Mathematical  Background  38 

3.  Derivation  of  Lattice  Filter  Equations-  43 

4.  Generalized  Analysis  (Lattice)  Filter--  56 

5.  Lattice  Filter  Advantages  62 

6.  Linear  Lattice  Filter  Applied  to 
Deconvolution  64 

III.   NONLINEAR  DECONVOLUTION  78 

A.  INTRODUCTION  TO  NONLINEAR  SYSTEMS  78 

B.  GENERALIZED  NONLINEAR  LATTICE  FILTER  81 

1.  Introduction 81 

2.  Nonlinear  Lattice  Filter  Development--  81 

3.  Nonlinear  Lattice  Applied  to 
Deconvolution  93 

4.  Inverse  Filtering  Simulation  Results--  96 
IV.   CONCLUSION 120 

APPENDIX  A:   LINEAR  LATTICE  FILTER  FORTRAN  PROGRAMS-  122 


APPENDIX  B:   NONLINEAR  LATTICE  FILTER  FORTRAN 

PROGRAMS 128 

LIST  OF  REFERENCES 139 

BIBLIOGRAPHY  141 

INITIAL  DISTRIBUTION  LIST  142 


ACKNOWLEDGEMENT 

I  wish  to  thank  Professor  S.  R.  Parker  for  his 
invaluable  support  in  the  writing  of  this  thesis.  His 
insight,  guidance,  and  encouragement  made  this  work  possible. 

My  thanks  also  go  to  Professor  L.  J.  Ziomek  for  his 
careful  review  and  editing  of  this  thesis.  His  efforts 
resulted  in  an  improved  final  product. 

Finally,  I  wish  to  thank  Major  Ed  Siomacco  for  sharing 
his  "lessons  learned"  in  thesis  preparation.  His  helpful 
hints  and  suggestions  were  greatly  appreciated. 


I.  INTRODUCTION 

T'he  problem  of  estimating  a  signal  based  upon 
observations  of  a  related  signal  is  one  of  the  most 
important  operations  in  signal  processing.  The  input  and 
output  signals  of  a  linear  system  are  related  by  the 
convolution  operation 


■/: 

-'  —  on 


y(t)  =  h(t)  *  x(t)  =J    h(t-T)x(T)dT-  (1.1) 

•ea 

where   h(t)  is  a  causal,  linear  time-invariant  (LTI),  system 

impulse   response.    Since  this  thesis  deals  primarily   with 

discrete   digital   signals,    continuous   time   signals   are 

sampled   at  uniform  intervals,   T,   and  are   represented   by 

discrete   sequences   x(nT)   =   x(n)   for   n   =   0,1,2,...,N. 

Discrete  convolution  for  a  LTI  system  is  defined  by 

n 
y(n)  =  h(n)  *  x(n)  =   £  h(n-ra)x(m)  (1.2) 

m=-oo 

and   is   shown   in  Figure  1.1.    The    basic    deconvolution 

problem  is  to  estimate   the  signal  x(n),  assuming  that   both 

y(n)   and   h(n)   are  known.    Figure  1.2  depicts  the  inverse 

filtering   process   of  recovering  the  input  signal  from   the 

output  signal  by  removing  the  system's  impulse  response. 

Deconvolution   has  important  applications  in  a   variety 

of  fields:    Radar,  communications,  image  processing,  speech 
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Figure   1.1:     Convolution 
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Figure  1.2      Deconvolution 


synthesis,  seismology.  For  example,  in  image  processing, 
deconvolution  is  used  to  recover  the  representation  of  the 
original  object,  x(m,n),  from  the  representation  of  its 
image,  y(m,n),  by  removing  the  blurring  caused  by  the  opti- 
cal system's  point-spread  function,  h(m,n).  A  common 
problem  which  arises  in  geophysics  involves  the  deconvolu- 
tion of  a  seismic  trace,  y(n),  into  the  approximately  known 
impulsive  waveform,  h(n),  and  the  desired  reflection 
response,  x(n),  which  reveals  the  structure  of  the  layered 
Earth.  In  other  applications,  h(n)  may  represent  the 
impulse  response  of  a  transmission  channel,  magnetic 
recording  medium,  or  measurement  device  which  broadens  and 
smears  ( intersymbol  interference)  the  desired  message  x(n). 

There  have  been  numerous  approaches  to  the  linear 
deconvolution  problem,  including  least-squares  filtering, 
linear  inverse  theory,  linear  programming,  and  homomorphic 
signal  processing.  The  first  part  of  this  thesis  deals  with 
least-squares  filtering  techniques;  the  application  of 
Kalman,  waveshaping,  and  lattice  filters  to  deconvolution  is 
reviewed.  Next,  the  lattice  filter  discussion  is  extended 
to  the  theory  of  the  generalized  lattice  filter.  Finally, 
nonlinear  system  theory  is  briefly  reviewed,  and  the 
nonlinear  lattice  filter  is  developed  and  applied  to  the 
inverse  filtering  problem.  Computer  simulation  results  for 
both  the  linear  and  nonlinear  lattice  filters  are  presented. 


II.  LINEAR  LEAST-SQUARES  DECONVOLUTION  TECHNIQUES 

A.   INTRODUCTION 

This  chapter  begins  with  a  discussion  of  the  principles 
of  least-squares  filtering  theory  since  deconvolution  is 
essentially  an  inverse  filtering  process.  Three  specific 
types  of  least-squares  filters  are  then  introduced:  Kalraan, 
waveshaping,  and  lattice  filters.  The  application  of  each 
of  these  three  filter  types  to  the  problem  of  linear  decon- 
volution is  studied.  Finally,  the  chapter  concludes  with 
the  presentation  of  computer  simulation  results  obtained  for 
deconvolution  experiments  using  the  lattice  filter. 

The  goal  of  deconvolution  is  to  recover  the  input 
signal,  x(n),  to  a  system  based  upon  observed  values  of  the 
system's  output  signal  y(n).  An  optimal  processor  must  be 
determined  to  produce  the  best  possible  estimate  of  x(n) 
based  upon  present  and  past  values  of  the  output  y(n).  A 
traditional  measure  for  defining  the  "best"  signal  processor 
is  the  minimum  mean-squared  error  criterion.  Using  this 
criterion,  the  estimate  x(n)  is  defined  by  a  linear 
combination  of  the  observed  values  y(n).  Assuming 
causality,  and  that  the  observed  signal  is  windowed  to 
include  only  the  M  past  samples,  x(n)  is  written  as 
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A  n 

x(n)  =  Y        h(n,i)y(i)    ,   n  >_  0  (2.1) 

i=n-M 

where  h(n,i)  =  0  for  n  <  i.  The  estimation  error  is  given  by 

e(n)  =  x(n)  -  x(n).    To  find  the  optimum  signal   processor, 

the   h(n,i)   coefficients   which   minimize   the   mean-square 

estimation  error  J,  where 

2  ^2 

J  =  E[e  (n)]  =  E[(x(n)  -  x(n))  ]   ,  (2.2) 

must   be   determined.    This  is  known  as  the   Wiener   filter 
formulation  of  the  problem.  [Ref.  1 : pp .  113,116] 

The  least-squares  theory  of  filtering  began  in  the 
1940's  with  the  work  of  Norbert  WIENER  [Ref.  2:pp.  147-148]. 
WIENER  developed  a  frequency  domain  procedure  to  design 
optimum  filters,  where  optimality  was  defined  by  minimizing 
a  mean-square  error  performance  criterion.  The  Wiener 
filter  is  conventionally  applied  to  linear  time-invariant 
systems  with  stationary  statistics  when  it  is  desired  to 
separate  one  signal  from  another.  In  the  early  1950's,  the 
Wiener  filter  was  extended  to  include  time  varying  and 
nonstationary  statistics,  but  the  calculations  are 
cumbersome  [Ref.  3:p.  1]. 

The  mean-square  estimation  error  J(n)  is  minimized  by 
setting  its  partial  derivatives,  with  respect  to  each  of  the 
filter  coefficients  h(n,i),   equal  to  zero.   Thus: 
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r)J(n) 

=   0  (2.3) 

0h(n,i) 


for  i  =  n-M,n-M+l , . . . , n.  This  yields  a  set  of  M+l  linear 
simultaneous  equations,  called  the  orthogonality  equations 
where 

R   (n,i)  =  E[e(n)y(i) ]  =  0  (2.4) 

ey 

for    n-M   <.   i   <_   n    and   n   >_   0.     R   (n,i)   is    the 

ey 
crosscorrelation   function  between  the  error   signal,   e(n), 

and  the  data  y(n).   Substituting  e(n)  =  x(n)  -  x(n)  into  the 

orthogonality   equations  gives  the  normal  ,   or   Wiener-Hopf 

equations : 

n 
E[x(n)y(i)]  =   Y   h ( n , k)E[y ( k ) y ( i ) ]  (2.5a) 

k  =  n-M 


or 


n 
R   (n,i)  =    £   h(n,k)R   (k,i)   .  (2.5b) 

xy         k=n-M        yy 


Since   the  autocorrelation  function  R    of  the  input   signal 

yy 

and   the  crosscorrelation  function   R   of  the  desired  output 

xy 
signal  with  the  input  signal  are  known  quantities,   the   M+l 

equations   can   be   solved  for  the   optimal   filter   weights 

h(n,i),  i=n-M,...,n.   If  data  vectors  are  defined  so  that 

T 
x(n)  =  [x(n-M),  x ( n-M+1 ),..., x ( n ) ]  (2.6a) 
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and 


y_(n)  =  [y(n-M),  y ( n-M+1 ),..., y( n ) ] 


(2.6b) 


then  the  M+l  Wiener-Hopf  equations  can  be  written  as 

T  T 

E[x(n)y  (n)]  =  hE[y_(n)y  (n)  ] 


(2.7) 


where  h  =  [h(n,n-M),  h( n-M+1 ),..., h ( n , n )] .  Assuming  that 
the  signals  x(n)  and  y(n)  are  stationary,  then  the  filter 
coefficients  are  time-invariant  and  h(n,k)  =  h(n-k);  and 
equation  (2.7)  can  be  written  in  matrix  form  as 


Ryy(O)  Ryy(l)  Ryy(2) 
Ryy(l)  Ryy(O)  Ryy(l) 
Ryy(2)   Ryy(l)   Ryy(O) 


Ryy(M) 


Ryy(M 


Ryy ( 0 ) 


h(0) 

—     , 

h(l) 

h(2) 

= 

• 

h(M) 

Rxy(O) 
Rxy(l) 
Rxy(2) 


Rxy(M) 


(2.8) 


Now,  the  optimal  coefficients  can  be  solved  by  inverting 
the  autocorrelation  matrix,  or  by  exploiting  the  matrix's 
Toeplitz  structure  (all  elements  are  the  same  on  any  given 
diagonal)  to  employ  the  more  efficient  Levinson  algorithm. 
(Levinson's  algorithm  will  be  discussed  in  the  development 
of  the  analysis  lattice  filter.) 

The   minimized  value  of  the  mean-square  estimation  error 
can  now  be  computed,  and  is  found  to  be 


2  n 

J    (n)  =  E[x  (n)]  -   Y     h(n,i)E[y(i)x(n) ] 

i  =  n-M 


(2.9) 


mm 
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T  T     -1 

=  R   (0)-E[x(n)y_  (n)  ]E[y_(n)y_  (n)]   E[£(n)x(n)] 
xx 

[Ref.   4:p.   148].    These  values  correspond  to  the  diagonal 

elements  of  the  covariance  matrix  of  the   estimation   error, 

T 
R    =  E[e(n)e  (n) ]  (2. 10) 

ee 

A 

where  e(n)  =  x(n)  -  x(n)  [Ref.  l:p.  120].  Furthermore,  the 
estimate  of  x(n)  is  given  by 

a  T  T     -1 

x(n)  =  E[x(n)y_  (n)]E[y(n)y_  (n)]   y_(n)  (2.11) 

=  hy(n) 

This  estimate  can  be  thought  of  as  the  projection  of  the 
desired  signal  x(n)  onto  the  space  spanned  by  the  components 
of  the  observation  vector  y(n).  The  minimized  estimation 
error  vector  is  orthogonal,  or  normal,  to  the  estimate 
x(n).  [Ref.  4:p.  147] 

This  completes  the  overview  of  least-squares  filtering. 
The  following  sections  present  several  linear  deconvolution 
techniques  which  employ  this  criterion:  Kalman  filtering, 
spiking  filters,  and  lattice  filters. 

B.   KALMAN  FILTER 

In  the  early  1960 's,  R.E.  KALMAN  introduced  an  optimal 
recursive  filter  based  on  state-space  time-domain  methods 
[Ref.  2:pp.  267-268].   The  Kalman  filter  estimates  the  state 
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of  a  linear  system,  and  is  optimal  in  the  sense  that  it  min- 
imizes the  mean-square  error  of  the  state  estimate.  The 
Kalman  filter  is  useful  when  the  system  is  defined  by  state 
space  equations:  The  system  signals  are  represented  by 
random  processes,  and  the  data  observations  are 
contaminated  by  noise.  The  Kalman  filter  algorithm  processes 
measurement  data,  and  requires  a  priori  state  space  models 
(known  or  assumed)  of  the  system  and  measurement  dynamics. 
Also,  the  statistics  of  the  system  input  and  measurement 
noises,  as  well  as  initial  condition  information,  are 
required  to  produce  the  state  estimate.  This  process  is 
depicted  in  Figure  2.1.  Here,  the  discrete  Kalman  filter 
equations  will  be  presented,  and  then  their  application  to 
deconvolution  will  be  described. 

The  state  space  representation  of  the  discrete  system 
and  measurement  models  (see  Figure  2.2)  are  written  as  [Ref. 
2:pp.  195-200] : 

x(k)  =  F(k-l)x(k-l)  +  w(k-l)  (2.12) 

z(k)  =  H(k)x(k)  +  v(k)  (2.13) 

where 

x(k)  =  (n  x  1)  system  state  vector 

F(k)  =  (n  x  n)  transition  matrix 

w(k)  =  (n  x  1)  system  noise  vector 

z(k)  =  (m  x  1)  measurement  vector 
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H(k)  =  (m  x  n)  observation  matrix 

v(k)  =  (m  x  1)  measurement  noise  vector. 

(Note  that  the  time  index  k  is  used  here  to  be  consistent 
with  the  discrete  Kalman  filter  literature.) 

The   noise  vectors  w(k)  and  v(k)  are  assumed  to  have   zero 

j 

J  mean,  white,  Gaussian  distributions  with  covariance  matrices 

.  T  T 

I  of    Q(k)    =    E[w(k)w  (k)]   and    R(k)    =    E[v(k)v  (k)], 

respectively.    Additionally,    w  and  v   are  uncorrelated  so 
]  T 

\  that  E[w(k)v  (j)]  =  0  for  all  k  and  . j  .   The  state  estimation 

J  error  is  defined  by 

> 

j 

•  e(k|k-l)  =  x(k)  -  x(k|k-l)  (2.14) 

and  the  associated  (n  x  n)  error  covariance  matrix  is 

T 
P(k|k-1)  =  E[e(k|k-l)e  (k|k-l)].  (2.15) 

An  updated,  a  posteriori  estimate  of  x(k)  is  obtained  from 
the  measurement  z(k)  and  the  a  priori  state  estimate 
x(kik-l)  by 

x(k)  =  x(k|k-l)  +  K(k) [ z ( k ) -H( k )x ( k I k-1 ) ]         (2.16) 

where  K(k)  is  the  (n  x  m)  Kalman  gain  matrix.  The  a 
posteriori  error  is  given  by 

e(k)  =  x(k)  -  x(k)  (2.17) 

and   the  associated  a  posteriori  error   covariance  matrix  is 
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T 
P(k)  =  E[e(k)e  (k)] .  (2.18) 


The   mean-square  estimation   error   criterion   is   minimized 

when 

T  T  -1 

K(k)=P(k|k-l)H  (k) [H(k)P(k|k-l)H  (k)+  R(k)]        (2.19) 

This  optimal  value  of  the  Kalman  gain  matrix  minimizes  the 
individual  terms  along  the  main  diagonal  of  P(k). 
Substituting  the  optimal  gain  matrix  K(k)  into  the 
expression  for  P(k)  results  in 

P(k)  =  (I_  -  K(k)H(k)  )  P(klk-l)  ,  (2.20) 

where  J_  is  the  identity  matrix.  In  order  to  compute 
equations  (2.16)  and  (2.19)  recursively,  the  a  priori 
estimates  x(k+l|k)  and  P(k+l|k)  must  be  determined  at  time 
k.   The  a  priori  estimates   are   given  by 

x(k+l|k)  =  F(k)x(k)  (2.21) 

T 
P(k+l|k)  =  E[e(k+l|k)e  (k+l|k)]  (2.22) 

T 
=  E[F(k)e(k)+w(k) ) (F(k)e(k)+w(k) )  ] 

T 
=  F(k)P(k)F  (k)   +  Q(k)  . 

The  Kalman  filter  algorithm  is  implemented  by  recursive- 
ly computing  equations  (2.16),  (2.19),  (2.20),  (2.21),  and 
(2.22).    Figure  2.3  is  a  diagram  of  the  Kalman  filter.   The 
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Figure  2.3      Discrete  Kalman  Filter 
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algorithm  is  initiated  with  the  initial  conditions 


x(0)  =  E[x(0) ]  (2.23a 


and, 

T 
P(0)  =  E[(x(0)-x(0) )(x(0)-x(0) )  ].  (2.23b) 

In  the  event  that  a  controlling  input  or  a  deterministic 
disturbance  u(k)  is  applied  to  the  system,  the  only  change 
in  the  above  algorithm  is  the  state  model  and  the  a  poster- 
iori estimate.   They  become  [Ref.  5:p.  130] 

x(k)  =  F(k-l)x(k-l)  +  w(k-l)  +  G(k-l)u(k-l)        (2.24) 

x(k)  =  F(k-l)x(k-l)+G(k-l)u(k-l)  (2.25) 

+K(k) [z(k)-H(k)F(k-l)x(k-l) ] . 

where  G(k)  is  a  (n  x  q)  matrix  and  u(k)  is  a  (q  x  1)  vector 
of  input  signals. 

The  problem  of  estimating  a  desired  signal  based  on 
noisy  data  observations  pertains  to  such  fields  as 
communications,  controls,  and  geophysics.  However,  the 
Kalman  filter  can  also  be  applied  to  these  deconvolution 
problems.  As  an  example,  the  discrete  Kalman  inverse  filter 
will  now  be  applied  to  an  exploration  seismology  problem. 

Typically,  in  the  search  for  underground  oil  and  gas 
deposits  ,  a  vibratory  signal  source  generates  a  pulse  of 
energy   which   is   transmitted  into  the   earth.    In   modern 
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seismology,  the  shape  of  this  source  wavelet  can  be 
carefully  controlled;  it  is  chosen  so  that  it  contains  only 
those  frequencies  which  are  transmitted  best  by  the  earth. 
As  the  source  wavelet  propagates  through  the  earth,  it 
encounters  many  different  layers  with  various  acoustic 
impedances.  At  these  layers,  both  partial  reflection  and 
refraction  occur,  creating  numerous  transmission  paths.  The 
received  seismic  signal  at  the  surface  is  composed  of  many 
overlapping  reflected  wavelets.  Therefore,  the  seismic 
trace  can  be  represented  as  the  convolution  of  the  original 
source  wavelet  with  an  impulse  train  representing  the 
various  layers  of  the  earth.  Moreover,  the  seismic  trace  is 
contaminated  by  measurement  noise  and  by  the  phenomena  of 
ghost  reflections  and  reverberations.  [Ref.  6:pp.  14-15] 
The  seismic  trace  can  be  described  mathematically  by 

z(t)  =  s(t,T)*r(t)  +  v(t)  =  /   s(t,T)r(T)dT  +  v(t)  (2.26) 

J    t0 

where   z(t)  =  measured  seismic  trace 

s(t,T)  =  finite  duration,  time  varying  wavelet 

r(t)  =  reflectivity  function  of  earth's  structure 

v(t)  =  measurement  noise. 

The  seismologist  must  extract  the  structure  of  the  earth, 
r(t),  by  analyzing  the  noisy  seismic  data  z(t).  This 
process  of  removing  the  wavelet  shape  from  the  trace  and 
leaving   behind  the  impulse  train  representing  the  reflected 
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wavelet's    strength    and    arrival    time    is    that    of 

deconvolution.    CRUMP  [Ref.   3:pp.   6-7]   shows  that  if  the 

seismic   signal   is  sampled  at   discrete,   uniformly   spaced 

intervals,   and  if  s(t,T)  and  r(t)  are  assumed  to  be  causal, 

then  z(t)  is  represented  by 

J 
z(k)  =  Y      [s(k,k-i+l )r(k-i+l > ]  +  v(k)  (2.27) 

i  =  l 

where  the  sample  number  k  =  1,2,3,...,  and 

L  =  length  of  the  wavelet  given  in  number  of  samples 

J  =  k  for  k<L 

=  L  for  k>_L  . 

When   M   traces   of  K  samples  in  length   are   available   for 

processing  then  z(k)  becomes  a  (M  x  1)  vector  where  the  j-th 

component  is  given  by 

J 
z  (k)  =   Y,    tH   (k)  )r(k-i  +  l  )  ]  +  v  (k)  (2.28) 

J        i=l   ji  J 

for  j  =  1,2, ... ,M  and  k  =  1,2,...,K.  This  assumes  that  the 
reflectivity  function  is  the  same  for  each  trace,  while  the 
shape  of  the  exciting  wavelet  may  vary  from  trace  to  trace. 
The  time-varying  wavelet  values  are  contained  in  the  (M  x  L) 
matrix  H(k)  where  the  j-th  row  contains  the  L  samples  of  the 
wavelet  which  generates  the  j-th  trace: 

H   (k)  =  s  (k,k-i+l )  .  (2.29) 

ji        J 

In  vector  form,  the  equations  become 
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z(k)  =  H(k)x(k)  +  v(k) 


(2.30) 


where  the  state,  measurement,  and  noise  vectors  are  given  by 

T 


x(k)  =  [r(k) ,r(k-l)  ,  .  . . ,r(k-L+l) ] 

T 
z(k)  =  [z  (k),z  (k) , . . . ,z  (k)] 
12  M 

T 
v(k)  =  [v  (k) ,v  (k)  ,  .  .  .  ,v  (k)]   . 
12  M  . 


(2.31a) 
(2.31b) 
(2.31c) 


This  is  the  Kalman  filter  measurement  model.    Now  the  state 
model  must  be  determined. 

The  state  model  is  arrived  at  by  assuming  a  general 
relationship  for  the  reflection  coefficients  of  x(k).  The 
assumed  relationship  is  the  autoregressive  equation 


r(k)  =  V      [b  (k-l)r(k-i)]  +  w(k-l)  . 
i  =  l    i 


(2.32) 


Comparing   equation  (2.32)  with  the  state  vector  x  ( k )  yields 
the  state  model 


x(k)  =  F(k,k-l)x(k-l)  +  w(k-l) 


(2.33) 


where   the  (M  x  L)  transition  matrix  and  the  (M  x  1)   system 
noise  vector  are  given  by 


F(k,k-1)  = 


b  (k-1)  b  (k-1) 
1        2 


b  (k-1) 
L 


0 


(2.34a) 
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T 
w(k-l)  =  [w(k-l) ,0,0, . . . ,0]   ,  (2.34b) 


X  is  the  identity  matrix,  and  0  is  the  null  vector. 
Equations  (2.30)  and  (2.32)  provide  the  measurement  and 
system  models  for  implementing  the  deconvolution  via  the 
Kalman  filter.  CRUMP  discusses  methods  by  which  to  obtain 
numerical  values  for  the  reflection  coefficient  vector  b(k) 
and  the  time-varying  wavelet  sample  matrix  H(k)  [Ref.  3 : pp 
8-11].  Once  these  matrices  are  determined,  the  recursive 
Kalman  filter  removes  the  effects  of  s(t,T)  from  r(t)  and 
generates  the  state  estimate  x(k)  which  provides  L  samples 
of  the  desired  reflectivity  function  at  each  time  k. 

It  is  interesting  to  note  that  both  the  Kalman  and 
Wiener  filters  are  minimum  mean-square  error  estimators, 
both  require  the  same  a  priori  knowledge  of  the  process  to 
be  estimated,  and  that  both  yield  identical  estimates. 
However,  the  Kalman  filter  does  have  distinct  advantages 
over  the  Wiener  filter.  First,  due  to  the  matrix  form  of 
the  state  space  equations,  the  Kalman  filter  has 
multichannel  capability  and  is  equivalent  to  a  bank  of 
optimal  estimators.  Moreover,  the  Kalman  filter  is  ideally 
suited  to  computer  implementation  due  to  its  discrete  and 
recursive  characteristics.  [Ref.2:pp.  268-269] 

The  discrete  least-squares  approach  which  follows  is  a 
viable    alternative    to   the   Kalman    filter    algorithm, 


25 


particularly   when   state   space   model   equations   are   not 
available  or  applicable. 

C.   LEAST-SQUARES  INVERSE  FILTER 

When  a  linear  system,   H(z),   is  excited  by  the  an  input 

signal  x(n) ,   the  output  y(n)  is  defined  by  the   convolution 

relationship   y(n)   =  x(n)  *  h(n).    Deconvolution   involves 

finding   the   inverse  filter  G(z)  such  that   H(z)G(z)   =   1. 

(Note   that  the  discrete  time  domain  is  related  to  the   fre- 

jwT 
quency  domain  through  the   equation  z  =  e     ,   where   T   is 

the  discrete   sampling  interval.)   This  condition  transforms 

to  h(n)  *  g(n)  =  d(n)  in  the  discrete  time  domain;   h(n)  and 

g(n)  are  the  impulse  responses  of  the  filters  H(z)  and  G(z), 

respectively,   and   d(n)  is  the  unit  impulse  function.     If 

this  condition  is  met,   then  the  original  input  signal   x(n) 

is  recovered  at  the  output  of  the  inverse  filter. 

The  transfer  function  of  the  causal  system  is  defined  by 

the   infinite   series   obtained  by  taking  the   one-sided   z- 

transform  of  the  system's  impulse  response: 

2       -i 
H(z)  =  Y(z)/X(z)  =  2L  Mi)z    .  (2.35) 

i  =  0 

H(z)   can   be  determined  by  inserting  a  known  sequence   x(n) 

into   the  system,   measuring  the  output  sequence   y(n),   and 

then   manipulating  their  z-transf orms .    The  inverse   filter 

G(z)   is   then   computed   by   carrying   out   the   polynomial 

division 


26 


-1       -2 
G(z)  =  1/H(z)  =  l/[h(0)+h( l)z   +h(2)z   +...]       (2.36) 

-1         -2  -M 

=  g(0)  +  g(l)z    +  g(2)z    +  ...  +  g(M)z    +  ... 


and   truncating   to  M+l  terms  if  necessary.    If   the   exact 

inverse   filter  is  approximated  by  truncating  G(z)  to   order 

M,   then  its  impulse  reponse  is  given  by  the  sequence   g(n), 

n  =  0 ,  1  ,  2 ,  .  .  .  ,  M.     Furthermore,    if   the   original   system's 

impulse  response  is  represented  by  h(n),  n=0 , 1 , 2 , . . . N,  then 

M 
d(n)  =  g(n)  *  h(n)  =  V  g(m)h(n-M)  (2.37) 

m  =  0 

for  0  <^  n  <_  N+M.   This  approximation  to  the  impulse  function 

improves  as  the  order  M  of  the  inverse  filter  is  increased. 

Now    the   stability   of   the   inverse   filter   will   be 

addressed.   If  H(z)  has  all  its  zeros  inside  the  unit  circle 

in  the  complex  z-plane,  it  is  referred  to  as  a  minimum-delay 

polynomial;   the   corresponding   sequence  h(n)  is   called   a 

minimum  phase-lag   sequence.   This  is  a  sufficient  condition 

to   guarantee  that  H(z)  has  a  stable   inverse,   because   the 

zeros  of  H(z)  become  the  poles  of  G(z),  and  G(z)  is  a  stable 

filter   if   all   its   poles  lie    within   the   unit   circle. 

Maximum-  and     mixed-delay    signals    are    obtained     by 

transforming   the  zeros  of  H(z)  from  z   to  1/z     where   the 

i        i 
superscript   '*'  represents  the  complex  conjugate  operation. 

A   maximum-delay  sequence  has  all  of  its  zeros   outside   the 
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unit  circle,  while  the  mixed-delay  sequence  has  zeros  inside 
and  outside  the  unit  circle. 

If   H(z)   has  N  zeros,   then   transforming   these   zeros 

N 
results  in  at  most  2   distinct  sequences.    Of  note,  each  of 

these   minimum,   maximum,   and  mixed-delay  signals  have   the 

2         * 
same  magnitude  spectrum:  |H(w)|   =  H(w)H  (w)  [Ref.  l:p.  98]. 

However,   they   do   have  distinct  phase  spectra   [Ref.   4:p. 

175].   The  maximum-delay  polynomial  can  be  written  as 

R        *        *       -1         *     -N 
H  (z)  =  h  (N)  +  h  (N-l)z    +...+  h  (0)z  (2.38) 

R 
The   so   called  reverse  polynomial  H  (z)   is   a   conjugated, 

reflected,   and   shifted  version  of  H(z).   The  corresponding 

R      *      *  * 

maximum  phase-lag  sequence  is  h   =  {h  (N),h  (N-l ),..., h  (0)} 

[Ref.  6:p.  72]   While  the  minimum-delay  filter  has  a  causal, 

stable   inverse  consisting  only  of  a   memory   function,   the 

maximum-delay  filter  has  an  inverse  which  consists  only  of  a 

stable,    noncausal ,    anticipation   function.    The   stable 

inverse   of   a  mixed-delay  function  consists  of  both   memory 

and   anticipation   functions.     Filters   with   nonvanishing 

anticipation  components  are  noncausal;   they  cannot  work   in 

real   time   since  the  future  values  of  the  filter  input   are 

not   available   for    processing.     This   problem   can    be 

circumvented  if  the  entire  signal  is  first  recorded  prior  to 

analysis;   then  the  required  future  input  data  is  available. 

[Ref.  4:p.  87] 
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The  energy  distribution  in  minimum,  mixed,  and  maximum 
delay  signals  will  now  be  examined.  Since  each  of  these 
signals  have  an  identical  magnitude  spectrum,  they  also 
have  the  same  total  energy.  However,  although  the  total 
energy  is  the  same,  the  rate  at  which  the  energy  builds  up 
differs  for  the  various  sequences.  Parseval's  theorem 
states  that  the  total  energy  in  a  signal  is  given  by 


L 


*  2  M        2 

|H(w)|   dw/(2n)  =  £  Ih(m)|   .  (2.39) 

-tr  m=0 

If  the  partial  energy  is  defined  as 

n        2 
P(n)  =  V    |h(m) I  ,  (2.40) 

m  =  0 

then   it  can  be  shown  that  the  energy  builds  up  quickest   in 

the   minimum-delay   sequence,   and   that  it   builds   up   the 

slowest  in  the  maximum-delay  sequence  [Ref.  6:pp.  75-76].  In 

other   words,   the  minimum-delay  signal  makes  its  impact   as 

soon   as   possible  since  its  energy  is  concentrated   at   the 

front   of  the  sequence.   The  maximum-delay  signal  makes   its 

major    impact    at   a   later   time   since   its   energy    is 

concentrated  at  the  end  of  the  sequence.    The  energy  curves 

associated   with   all  the  possible  mixed-delay   signals   lie 

between  these  two  extremes.    Finally,   it  can  be  shown  that 

the   convolution   of  two  minimum-delay   sequences   with   one 

another  results  in  a  minimum-delay  sequence.  The  convolution 

of    maximum-delay    signals   results   in   a    maximum-delay 
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sequence.  Convolution  involving  any  other  combination  of 
sequences  results  in  a  mixed-delay  sequence.  Of  note,  the 
resulting  maximum  and  minimum-delay  sequences  are  the 
reverse  of  each  other.  [Ref.  6 : pp .  73-74] 

The  method  described  above  of  finding  an  approximate, 
finite  length  inverse  filter  consisted  of  simply  truncating 
the  exact  inverse  filter  found  by  polynomial  division.  It 
was  seen  that  the  inverse  filter  G(z)  attempted  to  transform 
the  impulse  response  of  H(z)  into  a  unit  impulse  located  at 
the  origin.  This  can  be  thought  of  as  an  attempt  by  G(z)  to 
undo  the  blurring  effect  of  H(z)  (i.e.,  H(z)  "blurs"  the 
impulse  d(n)  into  the  impulse  response  h(n)).  If  the  input 
to  H(z)  is  designated  x(n),  and. if  the  output  of  G(z)  is 
x(n),  then  the  error  of  the  approximated  inverse  filter  is 

e(n)  =  x(n)  -  x(n)  for  0  <  n  <  N+M.  (2.41) 

The  error  energy  is  defined  by 

N+M   2 
J  =  V      e  (n)  (2.42) 

n  =  0 

For  the  polynomial  division  /  truncation  method,  J  decreases 

as  the  order  M  of  G(z)  increases  [Ref.  6:p.  136].  Seeking  to 

minimize   the   error  energy  J  with  respect   to   the   inverse 

filter  coefficients  leads  to  another  approach  for  finding  an 

approximate  inverse  filter. 


.M 
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Filters  which  minimize  the  mean-square  error  J  are 
called  least  error  energy  or  least-squares  inverses.  In 
general,  the  desired  output  sequence  x(n)  of  the  inverse 
filter  G(z)  could  be  of  any  shape.  G(z)  is  then  called  a 
waveshaping  filter.  When  applied  to  deconvolution,  the 
desired  output  of  the  inverse  filter  is  a  unit  impulse 
d(n-i),  where  i  =  0 , 1 , 2 ,  .  .  .  ,N+M  defines  the  lag  of  the  digi- 
tal inverse  filter.  In  this  case,  G(z)  is  called  an  i-th 
delay  spiking  filter  since  it  tries  to  condense  the  system 
impulse  response  h(n)  into  a  spike  with  i  delays.  Since  i  = 
0  ,  1 ,  2  ,  .  .  .  , N+M,  there  are  N+M+l  possible  spiking  filters.  As 
will  be  discussed,  there  are  preferred  values  of  the  lag  i 
for  minimizing  J,  depending  on  the  phase  characteristics  of 
h(n).   The  spiking  filter  problem  is  shown  in  Figure  2.4. 

It  is  convenient  to  restate  the  deconvolution  problem  in 
a  matrix  form  based  on  the  Yule-Walker,  or  autocorrelation, 
method  [Ref.  l:pp.  243-245].  This  is  accomplished  by 
defining  the  (M+N+l  x  M+l)  system  output  matrix  Y,  the 
(M+l  x  1)  impulse  response  vector  £,  and  the  (N+M+l  x  1) 
input  vector  x  as  follows: 
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(2.43a) 


g_  =  [g(0) ,  g(l) , 


,g(M)] 


(2.43b) 


x  =  [x(0) ,  x(l) ,  ...  ,x(N+M)] 


(2.43c) 


The  above  notation  is  for  the  general  case  of  the 
waveshaping  filter.  For  the  special  case  of  the  spiking 
filter,  .  the  y(n)  components  of  the  system  output  matrix  Y 
are  replaced  by  the  impulse  response  h(n).  Also,  the 
desired  signal  x  reduces  to  an  impulse  d(n-i).  Now 
equations  (2.37),  (2.41),  and  (2.42)  can  be  rewritten  as 


x  =  Yg. 


e  =  x  -  x  ,  and 


J  =  e  e 


(2.44a) 
(2.44b) 
(2.44c) 


As   discussed  in  the  introduction  to  this  chapter,   the 
least-squares   criterion  is  satisfied  by  minimizing   J   with 
respect  to  the  inverse  filter  coefficients  G.    This  results 
in  the  normal  equations 
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T       T 
Y  Yg_  =  Y  x  (2.45) 

T 
which   can   be   solved  for  g_.    By  recognizing  that   Y  Y   is 

equivalent  to  the  (M+l  x  M+l)  sampled  autocorrelation  matrix 

T 
R  =  E[y_£  ]  where  the  M+l  length   data  vector  is  y_  =  [   y(0), 

T  T 

y( 1 ) , . . . ,y(N) ,0,0, . . .0]  ,   and   that   Y  x  is  a  length   (M+l) 

cross-correlation  column  vector  r,  it  can  be  seen  that 

T   -1  T      -1 
£  =  (Y  Y)   Y  x  =  R   r  (2.46) 

-1 
where   R     can   be   evaluated    efficiently   by   Levinson's 

algorithm.  The  actual  filter  output  is  then 

-1  T 
x  =  Yg_=(YRY)x  =  Px  (2.47) 

T   -1  T 
where  P  =  Y(Y  Y)   Y   is  a  square  (N+M+l)  dimensioned   matrix 

called   the  projection  or  performance  matrix.    The  filter's 

performance   improves  as  P  approaches  the   identity   matrix. 

Now,  the  estimation  error  and  cost   function   are  written  as 

e  =  x  -  x  =  x(I  -  P)  (2.48) 

T      T 
J=ee=x(I-  P)x  (2.49) 


Since   in   the  deconvolution  problem  x  is  the  spike   with   i 

delays,   the   inverse  filter  output  x  is  actually   the   i-th 

-1  T 
column   of   the   P  matrix.    Also,   since   g_   =   R   Y  ,   the 

coefficients  of  the  i-th  spiking  filter  are  contained  in  the 

-1  T 
i-th   column   of  the  matrix   R   Y  .    Moreover,   the   energy 
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error  of  the  i-th  spiking  filter  reduces  to  J  =  1  -  P(i,i). 
Therefore,  in  order  to  realize  the  best  inverse  function 
(i.e.,  minimize  J),  select  the  delay  i  for  which  P(i,i)  is 
largest.  For  a  chosen  lag  i,  the  filter's  performance  also 
improves  as  the  order  M  of  the  inverse  filter  is  increased. 

Now  let  J(i)  represent  the  estimation  error  for  the  i-th 
spiking  filter.  The  grand  sum  of  squared  errors  is  defined 
as  V  =  J(0)  +  J(l)  +...+  J(M+N).  It  can  be  shown  that  V  = 
(M+N+l)  -  (M+l)  =  N,  where  N  is  the  order  of  the  system  H(z) 
[Ref.  4:p.  198].   Therefore,  V  is  independent  of  order  M. 

For  sufficiently  long  spiking  filters,  the  optimal  value 
of  the  delay  i  depends  on  the  phase  characteristics  of  the 
signal  h(n)  and  the  choice  of  the  lag  is  governed  by  the 
following  rules.  If  h(n)  is  a  minimum-delay  input,  the 
spiking  filter  should  have  zero  delay,  i  =  0.  This  says 
that  for  a  signal  with  its  energy  concentrated  towards  the 
front  of  the  sequence,  it  is  easiest  to  condense  it  to  a 
unit  impulse  at  the  origin.  If  the  signal  is  a  maximum- 
delay  input,  the  maximum-delay  spike  i  =  N+M+l  should  be 
selected.  If  h(n)  is  a  mixed-delay  signal,  then  the  i 
corresponding  to  the  largest  P(i,i)  should  be  chosen.  [Ref 
6.:p.  152] 

Under  certain  conditions,  the  estimation  error  will  go 
to  zero  as  the  length  of  the  inverse  filter  tends  to 
infinity.    First,   as   previously   discussed,   if  h(n)  is  a 
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minimum-delay  sequence  then  an  exact  inverse  can  be  found 
through  polynomial  division.  The  zero-delay  spiking  filter 
approaches  this  exact  inverse  filter  G(z)  as  the  number  of 
terms  M+l  increases.  Therefore,  the  estimation  error  goes 
to  zero  as  M  goes  to  infinity.  The  second  condition  is,  if 
h(n)  is  not  a  minimum-delay  sequence,  J  will  approach  zero 
as  1/M  if  the  lag  i  of  the  spiking  filter  is  chosen  to  be 
sufficiently  large.  [Ref.4:pp.  200-201] 

Up  until  now,  the  effects  of  measurement  noise  and 
imperfect  knowledge  of  the  distorting  function  H(z)  have 
been  neglected.  If  noise  is  introduced,  then  the  output 
y(n)  of  the  system  H(z)  driven  by  input  signal  x(n)  becomes 

y(n)  =  h(n)*x(n)  +  v(n)  (2.50) 

where  v(n)  is  assumed  to  be  zero-mean  white  noise  with 
variance  Q.  If  this  signal  is  passed  through  the  previously 
determined   inverse  filter  G(z),   the  filter  output   becomes 

x(n)  =  g(n)*y(n)  =  g ( n ) *h ( n ) *x ( n )  +  g(n)*v(n)      (2.51) 

=  d(n)*x(n)  +  u(n)  =  x(n)  +  u(n) 

where  u(n)  =  g(n)*v(n)  is  the  filtered  noise  signal  and  d(n) 

is   the   unit   impulse  function.    The  variance  of   u(n)  is: 

2  M   2         T 

E[u  (n)]  =  Q  £  g  (n)  =  Qg_  g.  (2.52) 

n  =  0 

[Ref.   l:p.  248].   This   variance   may   be  larger   than   the 
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original  noise  variance  Q.   For  example,   if  h(n)  is   a   low 

frequencey  signal,   then  g(n)  must  be  very  spiky  in  order  to 

compress   h(n)   into  a  high  frequency  content  spike.    As   a 

result,   g(n)   will   likely   have  values  greater   than   one. 

Therefore,   the  variance  of  u(n)  will  be  larger  than  Q.  This 

further  degrades  the  estimate  x(n). 

To  compensate  for  this  filtered  noise,   the  minimization 

criterion  is  modified  so  that 

N+M  2       T 

J  =   £   (d(n)  -  g(n)*h(n))   +  AQg.  g_  (2.53) 

n  =  0 

where  A  is  a  positive  parameter.   The  first  term  of  the  cost 

function  tries  to  produce  a  good  inverse  filter  whereas   the 

second  term  tries  to  reduce  the  output  noise.  If  A  is  large, 

noise   reduction  is  emphasized  at  the  expense  of  obtaining  a 

good  inverse  function.   A  small  A  emphasizes   finding  a  good 

inverse   function   g(n),   and  there  is  little   output   noise 

reduction . 

Using   this   new  cost   function,   the   resulting   normal 

equations  are  given  by 

T  T 

(Y  Y  +  AQI)g.  =  Y  x  .  (2.54) 

Comparing  this  with  equation  (2.45)  reveals  that  the  main 
diagonal  elements  of  the  sampled  autocorrelation  matrix  R 
have  been  modified  by  an  additive  term:  R(0)  becomes  R(0)  + 
AQ.    If  the  Backus-Gilbert  "prewhitening"  parameter  epsilon 
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is  introduced,  where 

£   =  AQ/R(0)  ,  (2.55) 

then  the  diagonal  elements  are  written  as  (1  +£  )R(0).   Even 

very   small   values  for  the  BG  parameter  have   a   beneficial 

T 
effect  of  stabilizing  the  inverse  of  the  matrix  (Y  Y  +  AQX) • 

[Ref.  l:p.  246] 

D.   LATTICE  FILTER 

1 .  Introduction 

The  lattice  filter  is  another  optimal  least-squares 
predictor  which  can  be  applied  to  the  linear  deconvolution 
problem.  The  lattice,  or  ladder  filter  derives  its  name 
from  the  cascade  form  of  its  signal  flow  graph.  Basically, 
the  lattice  filter  provides  an  alternative  to  the 
transversal  filter  for  modeling  a  signal.  It  can  be  viewed 
as  a  Gram-Schmidt  orthogonalization  of  the  incoming  data. 
This  will  be  elaborated  upon  in  subsequent  paragraphs.  To 
date,  much  of  the  work  done  with  lattice  filters  has  been  in 
the  areas  of  speech  analysis/synthesis,  seismology,  and  in 
high-resolution  spectral  estimation  [Ref.  7:p.  841].  Prior 
to  developing  the  lattice  filter  equations,  key  mathematical 
concepts   which  apply  to  this  discussion  will   be   reviewed. 

2 .  Mathematical  Background 

A  central  concept  behind  the  development  of  the 
lattice   filter  is  that  of  orthogonality.    Two  vectors   are 


38 


orthogonal  if  their  inner  product  is  zero.  The  geometric 
interpretation  of  this  is  that  the  vectors  are  at  right 
angles  to  one  another.  As  an  example,  if  the  random  varia- 
bles u  and  v  are  the  basis  vectors  of  a  two-dimensional 
space,  they  are  orthogonal  if  and  only  if  their  inner 
product  <u,v>  =  E[uv]  =  0.  In  the  case  where  the  linear 
space  is  spanned  by  the  random  variables  of  a  data  vector, 
two  vectors  are  orthogonal  if  the  corresponding  random 
variables  are  uncorrelated  and  if  one  or  both  have  zero 
mean  [Ref .  8:p.  92] . 

A  related  theorem  is  the  orthogonal  decomposition 
theorem,  which  states  that  any  random  variable  may  be 
decomposed  uniquely  with  respect  to  a  subspace  S  into  two 
mutually  orthogonal  parts,  one  part  which  is  parallel  to  S 
(i.e.,  lies  in  S)  and  the  other  part  orthogonal  to   S.   That 

,  A  A 

is,   y_   can  be  written  y_  =  y_  +  e  where  e  is  orthogonal  to   y_ 

and   to  the  basis  vectors   which  span  the   subspace   S,   and 

where   y_   is   defined  by  a  linear  combination  of   the   basis 

vectors.     If    S    is   spanned   by    the    basis    vectors 

M 

A       r— • 

{u( 1 ), u( 2 ),..., u( M) } ,   then  y  =   )   a(i)u(i)  where  it  can  be 

i  =  l 
shown   that   the   coefficients    are   given   by  the  equation 

2     -1 
a(i)  =   E[yu(i)]E[u  (i)]   .  [Ref.  l:p.  13] 

The   orthogonal  projection  theorem  adds  to   this   by 

stating:    the   orthogonal  projection  y_  of  a  vector  y  onto  a 

linear  subspace  S  is  that  vector  in  S  which  lies  closest   to 
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y_  with  respect  to  the  distance  induced  by  the  inner  product 
of  the  vector  space  [Ref.  l:p.  15].  Restated,  this  simply 
means  that  y_  is  the  best  estimate  of  y_  that  can  be  made  by  a 
linear  combination  of  the  basis  vectors  of  S  in  a  minimum 
mean-squared  error  sense.  Figure  2.5  shows  the  projection 
of  y_  into  the  subspace  S . 

Another  important  concept  in  developing  the  lattice 
filter  is  that  of  Gram-Schmidt  orthogonalization .  The  Gram- 
Schmidt  procedure  is  a  recursive  orthogonalization  process 
which  generates  a  set  of  mutually  orthogonal  basis  vectors 
{u ( 1 ), u ( 2 ),..., u ( M) }  from  a  given  set  of  basis  vectors 
{y ( 1 )» y( 2 ),..., y( M) } .  The  procedure  is  initialized  by 
letting  u(l)  =  y(l).  Next,  y(2)  is  decomposed  with  respect 
to  u(l).  That  part  of  y(2)  which  is  orthogonal  to  u(l) 
becomes  u(2).  Next,  y(3)  is  decomposed  with  respect  to  the 
subspace  spanned  by  (u(l),u(2)}.  That  part  of  y(3)  which  is 
orthogonal  to  this  subspace  becomes  u(3).  In  general,  the 
new  set  of  orthogonal  basis  vectors  is  defined  by 

n-1 
u(n)  =  y(n)  -  V   b(n,i)u(i)  (2.56) 

i  =  l 

2     -1 
for   2  <_  n  <_  M  and  where  b(n,i)  =   E  [  y  (  n  )u(  i  )  ]E  [u  (i)] 

Using    the   orthogonal   decomposition   theorem,    this    is 

equivalent   to   y(n)   =  y(n)  +  u(n)  where  y(n)  is   the   best 

estimate   of  the  n-th  component  of  the  data  vector  based   on 

the    previous   n-1   components   and   u(n)   represents    the 
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u 


Figure  2.5      Projection  of  y  onto  Space  S 
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estimation  error.   If  b(n,n)  =  1  is  introduced,  then 


n 
y(n)  =  V   b(n,i)u(i)   for  1  <  n  <_  M. 

This  can  be  written  in  matrix  form  as 


(2.57a) 


y_  -    Bu    where 


(2.57b) 


Z   =  [yd),  y(2),  ....  y(M)] 


u  =  [u(l) ,  u(2) , . . . ,  u(M)] 


B  = 


1  0 
b(2,  1)  1 
b(3,  1)    b(3,2) 


b(M,l)    b(M,2) 


This    is    a   convenient   notation   for   representing    the 

transformation  from  a  set   of  correlated  basis  vectors  y_   to 

an   uncorrelated  set  of  basis  vectors  u.   The  bases  y_  and   u 

span   the   same  M-dimensional  subspace  S,   but  there  are   no 

redundant   correlations  between  the  basis  vectors  of  set   u. 

Since  the  basis  of  u  is  uncorrelated,   its   components   u(i) 

(i  =  1,2,...,M)  are  referred  to  as  innovations  because   each 

additional  component  contributes  completely  new   information 

to  the  estimate  of  y_.  [Ref.  1 : pp .  16-18] 

Finally,   the  transformation  y_  =  Bu  corresponds  to  a 

LU   (lower  upper ) -Cholesky  factorization  of  the   correlation 

T 
matrix  of  y.    By  definition  R    =  E[yy  ].   Substituting  for 

yy 


42 


jr  results  in 

T   T         T 
R    =  BE[uu  ]B   =  BR   B  (2.58) 

yy  uu 

T 
where  B  is  lower  triangular,  B   upper  triangular,  and  R    is 

uu 
a    diagonal   matrix   since   the   basis   vectors   u(i)    are 

uncorrelated  with  one  another. 

3 .   Derivation  of  Lattice  Filter  Equations 

Now   the  lattice  filter  order  update  equations   will 

be   developed.    A  random  signal  y(n)  can  be  modeled  as   the 

output   of   a   causal,  stable,  linear  filter  H(z)   which   is 

driven  by  a  stationary,   uncorrelated,   white  noise  sequence 

{e( 1) ,e(2) , . . . ,e(P) }    [Ref.    l:p.    30].    A   P-th    order 

autoregressive  model  is  defined  by  H(z)  =  1/A  (z)  where 

P 
-1       -2  -P 

A  (z)  =  l+a(l)z   +a(2)z    +...+a(P)z     .      (2.59) 
P 

The  filter  A  (z)  has  several  names:   prediction  error  filter, 

P 
inverse  filter,  or  analysis  filter.   The  signal  y(n)   can  be 

written  as 

P 
y(n)  =  e(n)  -  T   a(i)y(n-i)   .  (2.60) 

1  =  1 

If  the  predictor  y(n)  is  introduced,   this  becomes 

y(n)  -  y(n)  =  e(n)  (2.61) 

P 
where  y(n)  =  -  >  a(i)y(n-i).   Now,   y(n)  is  the  estimate  of 

i  =  l 
y(n)  at  time  (n-1)  based  on  the  previous  P  samples,  {y(n-P), 
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y ( n-P+1  )  ,  .  .  .  ,  y(n-l)}  and  e(n)  is  equivalent  to  the  estima- 
tion error.  The  estimate  y(n)  can  be  thought  of  as  the 
projection   of  the  ( P+l ) -dimensional  y(n)  vector  onto  the  P- 

dimensional   subspace  spanned  by  the  components  of  the   data 

T 
vector   Y  =  [y ( n-1 ), y( n-2 ),..., y( n-P ) ]  .    In   general,   the 

past  values  of  y(n)  are  correlated  with  one  another.  There- 
fore, the  random  variable  components  of  Y  do  not  generally 
form  a  set  of  mutually  orthogonal  basis  vectors.  The  Gram- 
Schmidt  orthogonalization  procedure  removes  these 
correlations  and  transforms  Y  into  a  set  of  mutually 
orthogonal  basis  vectors  which  span  the  same  subspace. 
Additionally,   the  predictor  coefficients  a(i)  are   replaced 

by   the   reflection  coefficients  K  .    This  results   in   the 

i 
lattice  filter. 

In   order   to  obtain    the   optimal   estimator,   the 

predictor   coefficients   which   minimize   the    mean-squared 

2 
prediction  error  E[e  (n)]  must  be  determined.    This  problem 

was   discussed  in   the   introduction   to  this   chapter.   The 

resulting  matrix  form  of  the  normal  equations  is 


R(0)  R( 1 )  R( 2) 
R( 1 )  R(0)  R( 1 ) 
R(2)  R( 1 )  R(0) 


R(P) 


R(P) 


R(0 


1 
a(  1  ) 
a(2) 


a(P) 


Ep 

0 

0 


(2.62) 


where   R(k)   =   E [ y ( n+k ) y ( n ) ]  and  E   is  the   minimized  mean- 

P 
squared   prediction  error  for  the  P-th   order   filter.    The 
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sample  autocorrelation  matrix  components  are  calculated  from 
the  length  N  data  sequence  by 

N-l-k 
R(k)  =  (1/N)  V        y(n+k)y(n)  (2.63) 

n  =  0 

for  0  <_  k  <_  P   and  where  P  <_   N-l.  [Ref.  l:p.  150] 

Now   there   are  P+l  equations  and  P+l  unknown   model 

parameters  {a ( 1 ) , a( 2 ) , . . . , a( P ) ;E  }.  The  P+l  equations  can  be 

P 
solved   by   inverting   the   autocorrelation    matrix.    This 

3  2 

requires  0(P  )  operations  and  0(P  )  storage  locations.    The 

P+l   equations   can   be  solved  more   efficiently   by   taking 

advantage   of  the   matrix's   Toeplitz    structure   by   using 

Levinson's    algorithm.    Levinson's   algorithm  reduces   the 

2 
required  number  of  operations  and  storage  locations  to  CMP  ) 

and   CMP),    respectively   [Ref.  1 : pp .   150-151].    Being   a 

recursive    procedure,    Levinson's   algorithm   permits   the 

calculation  of  the  (P+l)-st  order  model  parameters  by   using 

the   previously  determined  P-th  order  model  parameters.   The 

matrix  form  of  the  algorithm  is  given  by 
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a    (1) 
P+l 

a    (2) 
P+l 


a    (P) 
P+l 

a    (P+l) 
P+l 


1 

a  (1) 
P 

a  (2) 
P 


a  (P) 
P 

0 


-   K 


P+l 


0 

a  (P) 
P 

a  (P-l) 
P 


a  (1) 
P 

1 


(2.64) 


where 


K 


£   a  (i)R(P+l-i) 
i  =  0   P 


P+l       P 

a  (i)R(i) 


i  =  0     P 


(2.65) 


The  P  subscript  on  the  "a"  parameters   specifies   the   P-th 

order   prediction   error   filter   A  (z)   whereas   the   index 

P 
identifies   the   appropriate  term  in  the   A  (z)   polynomial. 

P 
K      is   called   the  (P+l)-st  order   reflection   or   PARCOR 

P+l 
(partial  correlation)  coefficient..    The  PARCOR   coefficient 

K      represents  the  true,   or  direct,   correlation   between 

P+l 
y(n-P-l)   and   y(n)   with  the  effects   of   the   intermediate 

variables  (i.e.,  y ( n-P) , y(n-P+l ),..., y( n-1 ) )   removed.    The 

recursive  form  of  Levinson's  algorithm  is  written  as 


a    (m)=  a  (m)  -  K    a  (P+l-m)  for  l<_m<P,     (2.66a) 
P+l       P        P+l  P 


and 
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a    (P+l)  =  -K      for  m=P+l.  (2.66b) 

P+l  P+l 

This   shows  that  the  reflection  coefficient  for  each  stage  P 

is   equal  to  the  highest  coefficient  of  A  (z).    There  is   a 

P 
one-to-one   correspondence   between  the  PARCOR   coefficients 

and   the   coefficients   of  the  transfer  function  A  (z).   The 

P 
transfer   function   or,   equivalently ,   the   autocorrelation 

T 
matrix   R   =  E[y(n)y  (n)]  uniquely  determine  the   reflection 

coefficients  [Ref.   7:p.  829].   Taking   the   z-transform   of 

equation  (2.66)   results  in 


A    (z)  =  A  (z)  -  K    B  (z)  (2.67a) 

P+l        P        P+l  P 


where 


-1        -2  -P 

A  (z)  =l+a  (l)z   +a  (2)z   + . . . +a  (P)z    ,     (2.67b) 
P         P         P  P 


and 


-P     -1 
B  (z)  =  z   A  (z   )  ,  (2.67c) 

P  P 

-1  -2  -(P-l)   -P 

=  a  (P)+a  (P-l)z   +a  (P-2)z   + . . . +a  (l)z       +z 
P      P  P  P 

Due  to  the  effective  folding  about  the  axis  and  the  shifting 

to  the  right  of  the  sequence  A  (z),   the  polynomial  B  (z)  is 

P  P 

called   the   reverse  of  A  (z).    Taking  the  reverse  of   both 

P 
sides  of  equation  (2.67)  yields 

-1 
B    (z)=zB(z)-K    A(z)  (2.68) 

P+l  P        P+l  P 
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Combining  equations  (2.67)  and  (2.68)  into  matrix  form  gives 


A    (z) 
P+l 

B    (z) 
P+l 


-K 


P+l 


-1 


-K    z 
P+l 
-1 
z 


A  (z) 
P 

B  (z) 
P 


(2.69) 


The    forward   prediction   error    associated    with 
predicting   y(n)  from  the  previous  P  samples  {y ( n-P ) , . . . y ( n- 

A 

1)}  is  written  as  e  (n)  =  y(n)  -  y  (n)  where  the  subscript  P 

P  P 

represents   the   filter  order.   In  z-domain   notation,   this 

becomes 


E  (z)  =  A  (z)Y(z) 
P        P 


(2.70) 


Now,   the  backward  prediction  error  r  (n)  is  introduced.   It 

P 
is   defined   as   the   error   in   predicting   (   or   actually 

smoothing)  y(n-P)  from   the  future  P   samples   {y ( n-P+1 ,  .  .  .  , 

y(n)}.   It  is  written  as 


r  (n)  =  y(n-P)+a  ( 1 ) y( n-P+1 )+... +a  (P)y(n)    (2.71a) 
P  P  P 

or  in  z-domain  notation  as 


E  (z)  =  B  (z)Y(z) 
P        P 


(2.71b) 


The   optimal   backward  predictor  coefficients   minimize   the 

2 
mean-squared  smoothing  error   E[r  (n)].   Since   the   optimal 

P 
backward   and   forward   predictor   coefficients   are   mirror 
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2  2 

images  of   each  other,   then   E[r  (n)]  =  E[e  (n)].   That  is, 

the   backward  and  forward  prediction  error  vectors  have   the 

same  norms  [Ref.  8:p.  101]. 

As  will  be  shown,   at  each  instant  n,   the   backward 

prediction   errors  are  mutually  orthogonal   (i.e.,   uncorre- 

lated   if   assumed   to   be   zero   mean)   [Ref.    l:p.   170]. 

Therefore,    it    is    the    backward    prediction    errors, 

r  ( n ) , p=0 , 1 , . . . ,M,   where  M  is  the  filter  order,   which  form 

P 
the   new   set   of  basis  vectors.    To  demonstrate   that   the 

backward  prediction  errors  are  mutually  orthogonal,   let  M=3 

and   write   the  corresponding  equations  for  P  =   0,1,2,3   in 

matrix  form: 


r  (n) 
0 

r  (n) 

1 

r  (n) 
2 

r  (n) 
3 


0 


0     0 


a  (  1)   1      0     0 
1 

a  (2)  a  ( 1)   1     0 

2  2 

a  (3)  a  (2)  a  (1)  0 

3  3      3 


y(n) 
y(n-l) 
y(n-2) 
y(n-3) 


(2.72a) 


or 


r(n)  =  Lyjn) 


(2.72b) 


Note  that  the  first  column  of  L  contains  the  negatives  of 
all  the  reflection  coefficients.  Now  examine  the  covariance 
matrix  of  r ( n ) : 
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T  T 

R    =  E[r(n)r  (n)]  =  E[L.£(n)  (Lz(n)  )  ]  (2.73) 

rr 

T      T         T 
=  LE[y.(n)y.  (n)]L   =  LR   L   . 

yy 

Since   the  normal  equations  are  satisfied  within  the   matrix 

products   above,   R     reduces  to  a   diagonal   matrix   which 

rr 
verifies  that  the  components  of  r ( n )  are  uncorrelated.   That 

is 

T 
R    =  LR   L   =  D  =  diag{E  ,E  ,E  ,E  }  (2.74) 

rr      yy  0   12   3 

where  E   is  the  minimum  value  of  the  mean-squared  prediction 

P  2 

error  which  is  given  by  E[e    (n)  •]  .    It  can  be  shown   that 

2  P+l 

E     =  ( 1-K    )E   where   the   recursion  is   initialized  with 

P+l        P+l   P  2 

the  value  given  by  E   =   R   (0)  =   E[y  (n)]  [Ref.  8:p.  105]. 

0      yy 
Therefore,   the    prediction   error  decreases  by  a  factor  of 

2 
(1-K    )  from  one  lattice  stage  to  the  next.   Now,  since  the 

P+l 
elements   of   r(n)   are   uncorrelated,   equation   (2.72)   is 

equivalent  to  the  Gram-Schmidt  orthogonalization  of  the  data 

vector   yjn).    Also  note  that  rewriting  equation  (2.74)   as 

-1   -T 
R    =  L   DL    corresponds  to  a  LU-Cholesky  factorization   of 

yy 

the  covariance  matrix  of  y(n). 

To  complete  the  derivation  of  the  lattice  recursion 
equations,  multiply  both  sides  of  equation  (2.69)  by  Y(z)  to 
get 
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E    (z) 
P+l 

E    (z) 
P+l 


-1 


1 

-K         z 

E    (z) 

P+l 

P 

K 

-1 
z 

E    (z) 

P+l 

P 

(2.75) 


These   equations   are  transformed  into  -the  time   domain   to 
obtain  the  final  result: 


e    (n)  =  e  (n)  -  K    r  (n-1 ) 
P+l        P        P+l  P 

r    (n)  =  r  (n-1)  -  K    e  (n)  . 
P+l        P  P+l  P 


(2.76a) 
(2.76b) 


These  equations,  which  are  recursive  in  both  time  and  order, 
define  the  signal  flow  graph  of  the  analysis  lattice  filter, 
shown  in  Figure  2.6.  For  a  given  time  instant  n,  the  equa- 
tions are  evaluated  recursively  in  order.    The  inputs   into 

the   first   stage  of  the  lattice  filter  are  e  (n)  =  r  (n)   = 

0        0 
y(n).    The   successive   stages   of   the  filter  develop   the 

successively   higher   order  forward  and  backward   prediction 

errors.    The  output  from  the  final  stage  yields  the  desired 

M-th  order  forward  prediction  error,  while   all  lower  order 

prediction   errors   are   available  at   intermediate   stages. 

Since  the  backward   errors  are   generated  from  the  y_  data 

vector,   the   analysis   lattice  filter   actually   implements 

the  Gram-Schmidt  orthogonalization;  all   that  is  required  to 

implement    the    lattice    are    the    reflection    factors 

{K  ,K  ,  .  .  .  ,K  }  . 
1   2       M 
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Figure  2.6      Analysis  Lattice  Filter 
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Equation   (2.76)   can  be  manipulated  to   obtain   the 

lattice  form  of  the  synthesis  filter  H(z)  =  1/A  (z): 

M 

e  (n)  =  e    (n)  +  K    r  (n-1)  (2.77a) 

P        P+l        P+l  P 

r    (n)  =  r  (n-1)  -  K    e  (n)  .  (2.77b) 

P+l        P  P+l  P 

Figure  2.7  shows  the  corresponding  signal  flow  graph.    When 

the   input  to  the  synthesis  filter  is  the  forward  prediction 

error   sequence  e  (n),   the  output  is  the  original   sequence 

P 
y(n).    In   order  for  the  synthesis  filter  to  be  stable   and 

causal,   all   M   zeros  of  the  prediction-error  filter   A  (z) 

M 
must  lie  within  the  unit  circle  in  the  complex   z-plane.    A 

necessary  and   sufficient  condition  for  all  the  zeros  to   be 

inside  the  unit  circle  is  that  the  magnitude  of  each  of   the 

reflection   coefficients   {K  ,K  ,...,K  }  be  less   than   one. 

12  M 

[Ref.     l:pp.     168-169] 

The   reflection   factors   can  be   evaluated   by   several 

methods.     The    various   methods  arise   due   to   different 

definitions   of   the  optimality  criterion.    The    criterion 

used   here  of  minimizing  the  mean-squared  prediction   errors 

leads  to  what  MAKHOUL  calls  the  forward  and  backward  methods 

[Ref.  9].   They  are,  respectively: 

2 
K       =  E[e  (n)r  (n-1)]  /  E[r   (n-1)]        (2.78a) 
P+l,e       P     P  P 

2 
K       =  E[e  (n)r  (n-1)]  /  E[e   (n)]   .       (2.78b) 
P+l,r       P     P  P 
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Figure  2.7      Synthesis  Lattice  Filter 
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2  2 

Assuming   stationarity ,   and  since  E[r  (n)  ]  =  E[e  (n)  ]  as 

P  P 

previously    discussed,    then   the   forward   and    backward 

reflection  coefficients  are  equal.    That  is  K     =  K       = 

P+l     P+l,e 
K      .   The   Schwarz  inequality  implies  that  !  K  J  <_   1   for 

P+l,r  P 

each  P=1,2,...,M  [Ref.  8:p.  104].  An  alternate  technique  for 

calculating   the   reflection  coefficients  is   the   geometric 

mean   method.    It   was  introduced  by  ITAKURA  and   SAITO   in 

their   work   of   developing  a  digital  filter   structure   for 

time-domain   speech   analysis   [Ref.  10].   The  corresponding 

PARCOR  coefficient  is  given  by 

E[e  (n)r  (n-1)] 
P     P 

K     =  (2.79) 

P+l         2         2        1/2 
{E[e   (n)]E[r   (n-1)]} 
P         P 

These  PARCOR  coefficients  are  guarenteed  to  have   magnitudes 
less  than  one  [Ref.  7:p.  840]. 

A  third  method  was  used  by  BURG  in  the  maximum- 
entropy  method  of  spectral  estimation  [Ref.  11].  This  is 
the  harmonic-mean  method.  The  harmonic-mean  method  seeks   to 

minimize   the   sum  of  the  forward   and   backward   prediction 

2  2 

error  variances,   E[e   (n)]  +  E[r   (n-1)].    Minimizing  this 

P  P 

sum  with  respect  to  the  reflection  coefficients  results  in 
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2E[e  (n)r  (n-1)] 
P     P 

K     =  (2.80) 

P+l        2  2 

E[e   (n)  +  E[r   (n-1 ) ] 
P  p 

Again,  the  Schwarz  inequality  verifies  that  K   has  magnitude 

P 
o  less   than  one,   guaranteeing  that  the  synthesis   filter   is 

J 

£  causal  and  stable  [Ref.  l:p.  189], 

u 

4 .   The  Generalized  (Analysis)  Lattice  Filter 

z 

iaJ 

2  The   preceding   development  of  the  analysis   lattice 

i 

u  filter   equations   assumed   that  the   data   sequence   {y(n)} 

> 

D 

13  represented  a  time  sequence  with  stationary   statistics.   In 

this   section,   a   more  general   linear   prediction   problem 

UJ 

^  will  be   considered.     No   special   assumptions   are    made 

3  concerning   the  data.    The  data  values   need  not  be  delayed 
x 

2  versions  of  each  other;   they  need  not  even  represent  a  time 

sequence.  This  is  the  approach  taken  by  LENK  in  developing 
the  generalized  order  update  equations  [Ref.  12:p.  85].  The 
resulting  generalized  form  of  the  lattice  filter  makes  it 
suitable  for  multidimensional  and  nonlinear  signal 
processing  applications. 

Definitions  associated  with  the  normalized  form  of 
the  generalized  lattice  filter  will  now  be  introduced. 
First,  the  components  of  the  length  (M+l)  column  vector 
[y  ],  representing  a  single  realization  of  the  random 
process  Y,  are  designated  by  y  ,  A=  0,1,..., M.  The  forward 
prediction   error  in   estimating  the  (n+l)-st   element   from 
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the  previous  N  elements  of  the  data  vector  is  written  as 

N       M    N       x 
e     =   V  h  (n+1 )y  ( 2.81 ) 

n+1    A=0   A 

where  the   length  (M+l)   row   vector  of   coefficients  is 

given  by 


N  N       N  N 

[h  (n+l)]=[0, . . . ,0,-h     ,-h    , . . . , -h  , 1 , 0 , . . . , 0 ] . ( 2 . 82 ) 
*  n-N+1   n-N+2      n 

The   backwards   prediction  error  associated  with   predicting 

n-N 
y     from  the  next  N  elements  of  the  data  vector  is  given  by 

N      M  ^N       A 
r    =  Y      h  (n-N)y  (2.83) 

n-N   A=0   A 

where  the   associated   coefficients  are  given  by  the  vector 

^N  ^N       N  ^N 

[h  (n-N)  ]  =  [0,  .  .  .0,  1  ,-h     ,-h     ,  .  .  .  ,  -Ti  ,  0  ,  .  .  .  ,  0  ]  .  (  2  .  84  ) 
^  n-N+1   n-N+2       n 

The  norm  of  the  forward  prediction  error  is  defined  as 

N  N   2   1/2 

I |e   M  =  [E{(e   )  }]    .  (2.85) 

n+1  n+1 

The   norm   of  the  backward  prediction  error  is  defined  in   a 

similar  manner.    Now,   the  normalized  forward  and   backward 

N-th  order  prediction  errors  are  defined  by 

N  N  N  M         N  x 

e         =    e         /Me  It     =      £      a    (n+l)y  ,  (2.86a) 

n+1         n+1  n+1  ^=0       * 


and 


N  N  N  M         N 


X 


=    r         /    Mr  ||    =      Y.      b    (n-N)y  (2.86b) 


r  u      ~y 

n-N         n-N  n-N  ^=0 
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I  I 


where   the   normalized   prediction   error   coefficients   are 
defined  as 


N  N  N 

a  Jn+1)  =  h  (n+1)  /lie   I  I 
A  n+1 


X 


(2.87a) 


and 


N         ^N  N 

b  (n-N)  =  h  (n-N)  /  | |r   II 
A         A  n_N 


(2.87b) 


[Ref.  12:pp.  85-87] 

Using  the  normalized  form  of  the  generalized 
Levinson  algorithm,  LENK  demonstrated  that  equations  (2.87) 
could  be  updated  recursively  through  the  relation 


N+1 
a  .   (n+1  ) 

N+1 

b    (n-N) 
A 


N+1 


N 
a  (n+1) 

N 
b  (n-N) 


(2.88) 


where   the   partial   correlation   ( PARCOR )   coefficient    is 

(2.89) 


n 


N   __N 
K    =   E{e    r    } 
N+1       n+1  n-N 


and  where 


n 


0<K      )  = 

N+1 


n    2 
1  -  (K     ) 
N+1 


-  K 


n 


-  K 


n 
N+1 

1 


N+1 


(2.90) 


The  recursion  is  started  for  order  N=0  with  the  initial  pre- 
diction error  values  for  A=0,1,...,M   given  by   the  vectors 
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0  n+1 

[a  (n+1)]  =  [0.....0,  l/l|y    |»  ,0,...,0]    (2.91a) 

0  n 

[b  (n)]  =  [0,...,0,  1/lly   II,  0,...,0]  .     (2.91b) 

Multiplying  both  sides  of  equation  (2.90)  by  the  data  vector 
[y   ]  yields  the  desired  error  order  update  equations: 


N+1 
n+1 
N+1 
n-N 


n 


N+1 


N 

n+1 
N 
n-N 


(2.92) 


The  corresponding  signal  flow  graph  for  a  single  lattice 
filter  section  is  shown  in  Figure  2.8.  Figure  2.9  depicts  a 
third  order  generalized  lattice  filter.  [Ref.  1 2 : pp .  94-99] 

LENK  then  proved  that  the  reflection  coefficients 
(i.e.,  PARCOR  coefficients)  could  be  calculated  directly 
from  the  data  vector's  correlation  matrix  by  utilizing  the 
generalized  Schur  algorithm.  In  LENK's  notation,  the 
components  of  the  correlation  matrix  are  written  as  R 
E{y^  y   }.   Then  the  reflection  coefficients  are  given  by 


r\ 


n 


K 


N+1 


n- 

a 

N 

■N 

(n+1) 

n- 

P 

^  N 

-N 

(n-N) 

(2.93) 


where 
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N 
n+1 


N 
n-N 


Q-CKn  ]    ) 

N+1 


2.     -1/2 


N+1 
n+1 


a-[Kn  ]2> 

N+1 


-1/2 


N+1 
n-N 


a.  Single  Stage  of  Generalized  Lattice  Filter 
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b.  Alternate  Representation 


Figure  2.8      Generalized  Lattice  Filter 
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Figure  2.9    3— rd  Order  Generalized  Lattice  Filter 
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J 


IS 
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O 
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:j 
I 


and 


X  N    X  M   N      r 

CX         (n+l)  =  E{e    y  }  =   £a  (n+l)R 
N+l  n+l       p=0 


\  N    X      M   N       pA. 

/J    (n-N)  =  E{"    y  }  =   V  b  (n-N)R   . 
*N+1  n-N       n=0 


(2.94a) 


(2.94b) 


Equations   (2.94)   can   be   updated   through   the   recursion 
relation 


Of    <n+1> 
N+l 


fl         (n-N 
^N+l 


-N) 


n 
(K    ) 
N+l 


A 

01     (n+l) 

N 

P     (n-N) 
^  N 


(2.95) 


To   start   the   recursion  at  order  N=0,   the  values   of   the 
parameters  for  ,\=0,1,...,M  are  given  by  the  vectors 


n+l   -1    (n+l)0   (n+l)l 


(n+l)M 


\Oi  (n+l)]=||y    II    [R 

0 


,R 


,  .  .  .  ,  R 


X  n   -1    nO    nl        nM 

[//  (n)]  =  | |y   II    [R    ,R    , . . . ,R    ] 

^0 


]  (2.96a) 


(2.96b) 


[Ref.  12:pp.  100-101] 

5 .   Lattice  Filter  Advantages 

The   lattice  filter  has  several  advantages   compared 

to  the  direct,   or  transversal,  form  of  the  prediction  error 

filter    A  (z).     First   of   all,    due   to   the    built-in 

M 
orthogonalization  incorporated   into  the  lattice,  successive 

lattice   stages  are  decoupled.    Therefore,   the   reflection 
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coefficients  K  ,  P=1,2,...,M  are  independent  of  the  filter's 

P 
final  order.   The  M-th  order  least  squares  prediction  filter 

contains   all   the  prediction  error  filters  of  lower   order. 

Restated,  the  first  P  sections  of  the  M-th  order  filter  form 

the  P-th  order  prediction  filter.   Therefore,  lattice  stages 

may   be  added  or  subtracted  from  the  existing  lattice  filter 

without    having   to   recalculate   the   already    determined 

reflection  coefficients.   In  contrast,  when  the  order  of  the 

transversal   filter   is  changed,   all  the   A  (z)   prediction 

M 
error   filter   coefficients   must   be   recalculated.    As   a 

consequence,   to  specify  all  filter  orders  up  to  M   requires 

storage  of   M(M+l)/2   predictor  coefficients,  while  only    M 

reflection  coefficients  need  to  be  stored.  [Ref.  7:p.  830] 

Another  advantage  of  the  lattice  over  the 
transversal  filter  is  that  the  output  of  each  lattice  stage 
is  a  least-squares  prediction  error.  Therefore,  if  the 
desired  order  of  the  predictor  is  not  known  in  advance,  the 
output  of  the  various  stages  can  be  monitored  to  determine 
what  filter  order  is  adequate.  [Ref.  8:p.  106] 

Due  to  the  quantization  inherent  in  implementing 
digital  filters,  round  off  noise  is  introduced.-  Studies 
have  shown  that  the  performance  of  the  lattice  filter  is 
far  superior  to  that  of  transversl  filters  for  finite  word 
length  computations.  Furthermore,  the  filter  zeros  are  less 
sensitive     to    quantization   errors   in   the    reflection 
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coefficients  than  for  quantization  errors  in  the  transversal 
filter  coefficients.    Finally,   since  the  magnitude  of   the 
reflection  coefficients  is  less  than  one,  it  simplifies  the 
task   of   establishing  the  quantizer  overload   point.   [Ref. 
8:p.  106] 

Another  advantage  the  lattice  filter  holds  over  the 
transveral  filter  is  that  the  lattice  filter  is  minimum 
phase  (i.e.,  all  its  zeros  are  inside  the  unit  circle  so 
that  the  synthesis  model  is  stable)  if  and  only  if  the 
magnitude  of  the  reflection  coefficients  is  less  than  one. 
There  is  no  comparable  test  for  the  transversal  filter 
coefficients  to  determine  whether  the  synthesis  model  is 
stable.  [Ref.  8:p.  106] 

6 .   Linear  Lattice  Filter  Applied  to  Deconvolution 

The   transfer   function   of  the   lattice   filter   is 

determined   by   the  reflection  coefficients.   In   turn,   the 

values    of   these   reflection   coefficients   are    uniquely 

determined  by  the  prediction  error  filter  transfer   function 

A  (z),  or  equivalently  by  the  autocorrelation  sequence  R(k), 

M 
k  =  0,1,. ...M  [Ref.   7:p.  829].   Similarly,  equations  (2.70) 

and   (2.71b)  indicate  that  the   transfer   -functions  from  the 

input  y(n)  to   the  outputs   e  (n)  and   r  (n)   are  A  (z)   and 

M  M  M 

B  (z),  respectively.   Therefore,  the  analysis  lattice  filter 

M 
is  equivalent  to  the  whitening  filter  A  (z);   that  is,  it  is 

M 
the  inverse  of  the  system  model  H(z)  =  1/A  (z).   This  is  the 

M 
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basis   for   using   the   lattice   filter   in  a   deconvolution 

application. 

In   order  to  recover  the  input  signal  x(n)  from   the 

system's   output  signal  y(n),   it  is  necessary  to   determine 

the   inverse   of  the  system   transfer   function   H(z).    The 

initial  step  is  to  conduct  an  "identification  experiment"  to 

estimate   the  system's  parameters  (either  the  autoregressive 

filter    coefficients   or   the   lattice   filter    reflection 

coefficients).    To  accurately  estimate  the  parameters,   the 

chosen   experimental  input  signal  must  be  sufficiently   rich 

in   frequency   content,    or   more   formally,    persistently 

exciting  so  that  it  excites  all  the  modes  of  the  system.    A 

sequence   is  said  to  be  persistently  exciting  of  order  n   if 

its   (n  x  n)   autocorrelation   matrix   is  nonsingular   [Ref. 

13:pp.  70-71],   By  using  the  techniques  presented  in  section 

D.3  of  this  chapter,   the  resulting  output  sequence  y(n)  can 

be   processed  to  yield  the  desired  reflection   coefficients. 

When   the  analysis  lattice  filter  is  implemented  with   these 

coefficients  embedded  in  its  structure,   the  lattice  becomes 

the  inverse  of  A  (z)  =  1/H(z).    This  procedure  is   depicted 

M 
in  Figure  2.10. 

The  application  of  linear  lattice  filters  to  decon- 
volution was  simulated  using  the  FORTRAN  programs  LININV, 
ATOCOR,  LEVIN,  AND  LATICE.  These  programs  are  provided  in 
Appendix  A.    The   simulation  was   conducted   for   a   second 
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Figure  2.10      Deconvolution  Using  Lattice  Filter 
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order,  autoregressive ,  LTI  system  defined  by  the  equation 
y(n)  =  x(n)  -  (0.6y(n-l)  +  0.08y(n-2))  .  As  previously 
mentioned,  modeling  the  system  was  the  first  step  in  the 
deconvolution  process.  In  order  to  identify  the  "unknown" 
parameters,  the  system  was  excited  by  a  zero  mean,  unity 
variance,  Gaussian  white  noise  sequence.  The  output 
sequence  y(n)  was  processed  by  subroutine  ATOCOR  to 
calculate   the  components  of  the  autocorrelation  matrix  R 

yy 

Then  subroutine  LEVIN  implemented  Levinson's  algorithm  to 
evaluate  both  the  prediction  error  filter  coefficients 
and  the  associated  reflection  coefficients  from  the  autocor- 
relation matrix.  The  actual  output  of  subroutine  LEVIN  is 
the  lower  triangular  L  matrix  described  in  equation  (2.72); 
the  i-th  row  of  L  contains  the  i-th  order  prediction  error 
filter  coefficients  listed  in  reverse  order,  and  the  first 
column  contains  the  negative  of  the  reflection  coefficients. 
Once  the  reflection  coefficients  corresponding  to  1/H(z) 
were  calculated,  they  were  embedded  in  subroutine  LATICE 
which  implemented  the  analysis  lattice  filter  equations 
(2.76).  With  the  inverse  filter  of  H(z)  now  available 
(i.e.,  the  analysis  lattice  filter),  H(z)  was  driven  by  an 
"unknown"  signal  x(n).  The  output  of  H(z)  was  then  fed  into 
the  lattice  filter.  An  approximation  to  x(n)  was  recovered 
at  the  forward  error  output  signal  from  the  last  stage  of 
the  lattice. 
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Simulations  were  run  both  with  and  without  the 
presence  of  measurement  noise.  Table  2.1  displays  the 
resulting  L  matrices  for  one  simulation.  The  case  of  no 
measurement  noise  is  shown  in  Figures  2.11  through  2.14. 
Figure  2.11  is  a  plot  of  the  input  signal  x(n),  Figure  2.12 
depicts   the  system  output  y(n),   and  Figure  2.13  shows   the 

A 

lattice  filter  output  x(n).  As  can  be  seen,  the  results  of 
the  inverse  filtering  were  excellent--the  x  and  x  curves  are 
identical.  This  is  verified  by  Figure  2.14  which  shows  that 
the  mean-square  error  between  x(n)  and  x(n)  is  nearly  zero. 
The  running  average  mean-square  error  was  calculated  using 
the  equation 


n  2 

MSE(n)  =  ~\J   (1/n)  £  (x(i)  -  x(i))    .  (2.97) 

The  simulation  was  then  repeated  for  a  nonzero 
measurement  noise.  Here,  the  added  measurement  noise  was  a 
zero  mean,  0.0025  variance,  Gaussian  white  noise  sequence. 
Using  the  same  input  as  for  the  previous  case,  the  outputs 
of  the  system  and  lattice  filter  are  shown  in  Figures  2.15 
and  2.16,  respectively.  Figure  2.17  is  a  plot  of  the  mean- 
square  error  between  x(n)  and  x(n).  The  results  in  Table  2.1 
show  that  even  with  the  presence  of  measurement  noise,  the 
system  parameters  were  still  accurately  identified.  The 
lattice  filter  recovered  the  basic  shape  of  x(n),  however, 
it    could   not    remove    the   additive   measurement   noise. 
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Therefore,  the  output  of  the  lattice  filter  is  a  noisy 
or  "fuzzy"  version  of  x(n).   Additional  filtering  is  required 

A 

to  remove  the  noise  component  of  x(n).  The  results 
presented  here  are  representative  of  those  obtained  for 
simulations  involving  other  stable,  autoregressive ,  LTI 
systems . 
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TABLE  2 . 1 
RESULTS  OF  MODELING  A  LINEAR  SYSTEM 

Modeling  Problem: 

-1        -2 
System:   H(z)  =  1/[1  +  0.6z    +  .08z    ] 

System  input:  Gaussian  white  noise,  N(0,1) 

Number  of  white  noise  realizations  used:   25 

Number  of  points  per  realization:   5,000 

Modeling  Results: 

a.   No  Measurement  Noise 


L  = 

1 
0.555 
0.077 

0 

0 
1.0 
.598 

0 
0 

1 

Order, P    Ep 

Kp_ 

0 

1 
2 

1.445 
0.999 
0.993 

-0.555 
-0.077 

b.   Measurement  Noise  is  N(0, 0.0025) 


L  = 

1 
0.555 
0.077 

0 

0 
1.0 
.598 

0 
0 

1 

Order 

l£   l£ 

Kp 

0 

1 
2 

1  .449 
1  .002 
0.996 

-0.555 
-0.077 

70 


oo* t  sro     os*o  sz'o     oco  s<ro-  os'o-  sz*o-  oo*x- 

(N)X 


Figure  2.11    System  Input,  x(n) 
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Figure  2.12    System  Output,  y(n) 
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Figure  2.13    Lattice  Filter  Output,  x(n) 
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Figure  2.14    Mean-Square  Error  Between  x(n)  and  x(n) 
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Figure  2.15    System  Output  Plus  Measurement  Noise 
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Figure  2.16    Lattice  Filter  Output,  x(n) 
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Figure    2.17         Mean-Square   Error   Between   x(n)    and   x(n) 
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III.   NONLINEAR  DECONVQLUTION 

A.   INTRODUCTION  TO  NONLINEAR  SYSTEMS 

While   the   previous  chapter  dealt   solely   with   linear 

l 
r 

systems   or   plants,   this  chapter  will  address  the   inverse 
,'f  filtering  problem  involving  nonlinear  systems.    The  chapter 

'i  starts   with   a   brief  introduction   to   modeling   nonlinear 

jii; 

systems.  Then  it  quickly  proceeds  to  extend  the  generalized 
linear  lattice  filter  results  of  section  II. D. 4  to  a 
nonlinear  analysis  lattice  filter.  The  nonlinear  lattice 
filter  is  discussed  in  detail.  To  conclude,  results  from 
numerous  simulations  involving  the  lattice  in  deconvolution 
applications  are  presented. 

System  linearity  is  defined  in  terms  of  the  principle  of 
superposition.  If  the  rule  by  which  the  system  transforms 
the  input  x(k)  into  the  output  y(k)  is  represented  by  the 
operator  T,  then  the  system  is  said  to  be  linear  if  and  only 
if 

T[ax  (k)  +bx  (k)]  =  aT[x  (k)]  +bT[x  (k)]  (3.1) 

12  1  2 

=  ay  (k)  +  by  (k) 
1         2 

for  arbitrary   constants  a  and  b.    The  system  is   nonlinear 

if   equation   (3.1)   is  not  satisfied.    Estimation   of   the 

parameters  of  a  nonlinear  system  is  a  complex  problem. 
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One   approach   to  the  parameter  estimation   problem   for 

nonlinear  filters  is  the  Bayesian  approach.  This  approach 
leads   to   various  approximate  solutions  for   the   parameter 

estimates.    In  the  Bayesian  approach,   the   parameters   are 

considered   to  be  random  variables.    The  parameter  vector  h 

is   considered   a  random  vector  with  a   probability   density 

function  p(h).   Introducing  the  observation   vector  y_(t)  and 

the   input   vector  x(t),   both  of  which  contain  data   up   to 

time  t,   the  a  posteriori  probability  density  function  for  h 

is  p(  h  |y_(  t )  ,x(  t )  )  .    One  possible  choice  for  the  estimate  of 

A 

the  h  vector  is  to  select  the  conditional  mean  h ( t )  = 
E[h  |  y_(  t )  j_x(  t )  ]  .  Selecting  this  as  the  estimate  minimizes  the 
variance    of   the   parameter   estimation   error.     Another 

A 

possible  choice  is  to  select  the  h(t)  which  maximizes 
p(  h  |  y_(  t )  ,x(  t )  )  .  This  most  likely  value  is  known  as  the 
maximum  a  posteriori  (MAP)  estimate.  Finally,  the  Bayesian 
approach  also  leads  to  the  extended  Kalman  filter  for 
nonlinear  state  estimation  problems.  [Ref.  13:pp.  32-41] 

Another  popular  method  for  modeling  nonlinear  systems  is 
based  on  the  Volterra  series.  The  series  is  named  after  the 
mathematician  Vito  VOLTERRA.  The  first  person  to  apply  the 
series  to  nonlinear  systems  was  Norbert  WIENER.  If  the 
"black  box"  approach  is  taken  towards  a  nonlinear,  time- 
invariant  system,  the  relationship  between  the  input  x(t) 
and   the   output   y(t)  can  be  represented   by   the   Volterra 
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series : 


■oo 
h    (T    )x(t-T    )dT 


rtja 
=  h 


■oo 

h    (T    ,T    )x(t-T    )x(t-T    )dT    dT 
,_co      2       1       2  1  2         12 

•oo 

h    (T    ,T    ,T    )x(t-T    )x(t-T    )x(t-T    )dT   dT   dT 
,_oo       3       12      3  1  2  3         12      3 


+  /     .  .  ./     h    (T 


+/     .../     h    (T    ,...,T    )x(t-T    )...x(t,T    )dT    . . . dT 
-co      '-co  J  n  1  n         1  n 

+       . . .  (3.2) 

where  n  =  1,2,...   and   h(T   ,...,T   )=0   for  any  T   <  0, 

n   1        n  j 

j   =   l,2,...,n  .   The  functions  h  (T  ,...,T   )   are   called 

n   1        n 
Volterra    kernels   of   the   system.     This   equation  is   a 

functional  series.   That  is,  it  performs  an  operation  on  the 

function  x(t)   which   results  in  a  number  for  y(t).   If  the 

n-th  order  Volterra  operator,  H  ,  is  introduced  where 

n 

y    (t)    =    H    [x(t)]  (3.3) 

n  n 

rco        roo 

=j       . . J     h    (T    ,...,T    )x(t-T    )...x(t-T    )    dT    ...dT    , 

J_/nl  n  1  nl  n 

-co       -'-oo 

then  the  series  can  be  expressed  in  operator  notation  as 


y(t)  =  H  [x(t)]  +  H  [x(t)]  +  ...  +  H  [x(t)]  +...    (3.4) 
1  2  n 

00  CO 

=    £    h  [x(tn  =   £    y  <t)  . 

n=  1   n  n= 1   n 
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Figure  3.1   is  a   graphic   representaion   of  this   equation. 

[Ref.   14:pp.  7-9]   The  H   operator   is  a   linear   operator, 

1 
H   is  a  quadratic  operator;  and,  in  general,  the  H   operator 

2  n 

involves  the  term  x(t)  to  the  n-th  power. 


B.   GENERALIZED  NONLINEAR  LATTICE  FILTER 

1 .  Introduction 

In  this  section,  it  will  be  demonstrated  how  the 
generalized  lattice  filter  introduced  in  section  II. D. 4  can 
be  applied  to  modeling  nonlinear  systems.  The  development 
will  start  by  looking  at  the  discrete  form  of  the  Volterra 
series.  Based  upon  this  series  an  alternate  tensor  notation 
representation  will  be  introduced.  The  results  of  section 
II. D. 4  will  then  be  extended  to  handle  a  two-dimensional 
field  of  data  which  will  result  in  the  generalized  nonlinear 
lattice  filter. 

2 .  Nonlinear  Lattice  Filter  Development 

The  discrete  form  of  the  Volterra  series  of  equation 
(3.2)  is  given  by 


y(k)  =  h   +    h  (n  )x(k-n  )  (3.5) 

0       11       1 

+   h  (n  ,n  )x(k-n  )x(k-n  )  +  ... 
2   12       1       2 


Using  LENK's  tensor  notation,  an  equivalent  form  of  equation 
(3.5)  is  given  by 
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Figure  3.1      Volterra  Series  Representation 
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y(k)  =  H  +  H  x   +H   xx   +  ...  ,  where        (3.6a) 

X'  T 

x  \k)  =  tx(k) ,x(k-l ) , . . . ,x(k-A;) x(k-N)]   (3.6b) 

for  A*,  =  0,1,...,  N.  As  an  example,  if  N=l  and  equation 
(3.6a)  is  truncated  to  the  three  terms  shown  above,  then 
this  equation  can  be  rewritten  as 

y(k)  =  H  +  H  x(k)  +  H  x(k-l)  +  H   x(k)x(k) 
0         1  00 

+H   x(k)x(k-l)  +  H   x(k-l)x(k) 
01  10 

+  H   x(k-l)x(k-l )  .  (3.7 ) 

11 

Based  on  this  tensor  notation,   LENK  introduced  an  alternate 

representation   for   the   nonlinear   system.     Instead    of 

defining   the  components  of  the  vector  x  (k)  as  the   present 

and   past   values  of  the  signal  x(k)  as  in  equation   (3.6b), 

the  components  are  redefined  in  terms  of  x(k)  raised   to  the 

Jy  power  for  A;  =  0,1,2,...,N.   That  is 

■v.  0      12  NT 

x(k)  =  [x  (k),x  (k),x  (k),...,x  (k)]  (3.8) 

2  NT 

=  [  1    ,x(k)  ,x  (k) , . . . ,x  (k)] 

for  ^*,  =  0,1,2,  ...  ,N.  Now  the  output  of  the  nonlinear  system 
is  defined  by 

y(k)  =  x  (k)...x  (k-N)H   ...  (3.9) 

Ao       Ah 
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for  ^  , A, , . . .  =  0,1,2,..  .,P.   The  variable  P  gives  the  order 

of  the  filter  and  N  defines  the  filter's  finite  memory.   The 

term   H  v  .  .  ...   plays  a  role  similar  to  that  of  the   Volterra 
Ao   Ah 

kernel;  it  can  be  considered  a  (P+l)-order  tensor.  To 
clarifly  the  meaning  of  this  notation,  the  following  example 
with  N=l  and  P=2  is  provided: 

y(k)    =    H.    ,k  °(k)x    '(k-1)    . 

*o  Ai  2 

=    H      +H      x(k)+H      x(k-l)+H      x(k)x(k-l)+H      x    (k) 
00      10  01  11  20 

2  2  2 

+    H      x    (k-l)+H      x    (k)x(k-l)+H      x(k)x    (k-1) 
02  21  12 

'!!  2  2 

j>   '  +    H      x    (k)x    (k-1  )     .  (3.  10) 

'm  22 

[Ref.    12:pp.    39-48] 

The  above  nonlinear  model  for  y(k)  is  a  moving 
average  (MA)  model:  the  output  is  defined  in  terms  of  the 
input  signal.  The  next  step  is  to  establish  an 
autoregressive  ( AR )  model  where  y(k)  is  estimated  in  terms 
of  its  past  values.  This  type  of  model  is  useful  when  the 
input  signal  is  not  readily  available.  An  autoregressive 
model  is  obtained  by  redefining  the  observation  vector  in 
terms  of  y(k)  vice  x(k).   Then  the  AR  model  is  given  by 

y(k)  =  y  (k-D  .  .  -y  N(k-N)H   .  .  ..  (3.11) 

A,   An 

for  A\  =  1,2, ...,P.  If  the  system  is  driven  by  a  white 
noise   signal   u(k),   and   if  the   system's   parameters   are 
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approximated  by  H  ^  ...   ,  then  the  estimation  error  is 

Ao         A,>i 

e(k)=y(k)-y(k)=[y  (k-1 ) . . . y  1 k-N )H   ...   +  u(k)] 

Ai    AN 

-  [  y  \k-l)...y  "k-NJH*  .  .  .A  ]   (3.12) 

Ai    ^n 

If  the  system  parameters  are  known  exactly,  then  the  estima- 
tion error  and  the  white  noise  sequences  are  equal--that  is 
e(k)=u(k).  One  note  concerning  the  nonlinear  AR  model: 
unlike  the  linear  AR  model  whose  stability  is  easily 
determined  (i.e.,  the  system  is  stable  if  all  its  poles  lie 
within  the  unit  circle  in  the  complex  z-plane),  it  is 
difficult  to  judge  for  which  class  of  inputs  the  AR 
nonlinear  system's  output  will  remain  bounded.  This  is 
because  the  order  of  the  nonlinearity  increases  with  time. 
[Ref.  12:pp.  68-69] 

Using  this  alternate  tensor  form  of  the  AR  nonlinear 
system,  along  with  the  generalized  lattice  of  section 
II. D. 4,  the  nonlinear  lattice  filter  will  be  developed.  To 
simplify  the  discussion,  the  filter's  finite  memory  will  be 
restricted  to  N  =  2.  Then  the  product  Y(k)  =  [y  '( k-1  )  y  \  k-2 )  ] 
forAv,A^  =0,1,..., P  forms  a  second  order  tensor,  or  a  two- 
dimensional  data  field  given  by: 


85 


1  III 


Y(k)  = 


y(k-2) 


. . .   y  (k-2) 


y(k-l)   y(k-l)y(k-2)  ...  y(k-l)y  (k-2) 

2         2  2       P 

y  (k-1)   y  (k-l)y(k-2). . .  y  (k-l)y  (k-2) 


P         P  P       P 

y  (k-1)   y  (k-l)y(k-2) . . .  y  (k-l)y  (k-2) 


(3.13) 


J 

::/c 


ig 


The  types  of  systems  that  this  nonlinear  lattice 
structure  can  model  exactly  are  of  the  form  shown  in  Figure 
3.2  .  The  nonlinear  combinations  block  forms  a  weighted  sum 
of  the  cross-products  and  powers  of  the  input  variables 
y(k-i),  for  i=l,2,...,N.  As  an  example,  with  N=2  and  P=2  the 
general  equation  for  the  system's  output  when  excited  by  the 
input  x(k)  is  given  by: 


y(k)  =  x(k)  -  {  H    +  H   y(k-l)  +  H   y(k-2) 

11     21  12 

2  2 

+  H   y(k-l)y(k-2)  +  H   y  (k-1)  +  H   y  (k-l)y(k-2) 
22  31  32 

2  2 

+  H   y  (k-2)  +  H   y(k-l)y  (k-2) 
13  23 

2       2 
+  H   y  (k-l)y  (k-2)  }  (3.  14) 

33 


It   is  also  assumed   that  the   system  is  time-invariant; 

otherwise  the  model  changes  with  time  k.   In  order  to  define 

2 
the   forward   and   backward  prediction   errors,   the   (P+l) 

elements  of  the  two-dimensional  data  field  must  by  converted 
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Figure  3.2      Nonlinear  System  Model 
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into  a  one-dimensional   sequence.    The  ordering   chosen   by 
LENK  is  listed  below: 

order  =  {( P, P ),( P-l , P ),( P, P-l ),...,( 0 , P) , (P, 0 ) , 
(P-1,P-1) , (P-2,P-1), (P-l,P-2), . . .  , (0,P-1), 
(P-1,0),  .  .  .,(1,1), (0,1), (1,0), (0,0)}  ,   (3.15) 

where  the  components  of  the  data  field  are  identified  by  the 

indices   (m,n)   for  m,n  =  0,1,..., P.    The  elements   of   the 

2 
ordered   set  are  numbered  consecutively  from  0  to   (P+l)  -1. 

The   notation  (m,n)-q  is  used  to  identify  the   q-th   element 

prior   to   the   element  (m,n)  as  referenced  to   the   ordered 

sequence.   This   notation   will   be   further   abbreviated  to 

mn-q.   Now   the   (q-l)-order,   normalized,    forward    error 

associated    with    predicting    the  value   of   the   element 

m       n 
y  (k-l)y  (k-2)  from  the  preceding  (q-1)  elements  is  given  by 

_q-i         q-i  A,     Az 

e     =  a    (m,n)y  (k-l)y  (k-2)  (3.16) 

mn      A»  Ax 

q-1 

for  Ai  ,Ai     =0  ,  1  ,  2  ,  .  .  .  ,  P  .   Note  that  the  coefficients  a  ->  ->  can 

A, At 

be   thought  of  as  the  components  of  a  second   order   tensor. 

q-1 
The   coefficient   a.  -v  is   equal  to   zero  when   the   indices 

A»Ajl 

{X[   »^2,)  do   not   correspond   to  the  (q-1)  elements  preceding 

(m,n)  in  the   ordered   sequence   (i.e.j  when  (  A|»Aa )  >  (m,n) 

or  when   (AV,A^.)  <.  mn-q  ).   Also,  when  (  A,  »A^ )  =  (m,n),  then 

q-1  q-1 

a    (m,n)  =  1  /  | |  e      | |   .  (3.17) 

mn  mn 
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Similarly,   the   normalized  backward  prediction  error   in 
estimating  y(mn-q)  from  the  next  (q-1)  points  is  given  by 


Jl"1  q-1         A,       Xi 

r     =  b     (mn-q)y   (k-l)y   <k-2) 
mn-q    AjAa. 


(3.18) 


q-1 


for  A|  ,  Ax   =0 » 1 » 2 , . . . ,P.    The  coefficients  b.  ,  (mn-q)   equal 

AiAa. 

zero  when  the  indices  (Aj,X^)  do  not  correspond  to  the  (q-1) 
elements  following  (m,n)  (that  is,  when  (A^A^)  <  (mn-q)  or 
(  A,.  A^)  >  (m,n)  ).   When  (Ai>A^)  =  (m,n),  then 


q-1 
b    (mn-q)  = 
mn-q 


q-1 
1/  Mr    || 

mn-q 


(3.19) 


[Ref.  12:pp.  145-146] 

Using  the  normalized  nonlinear  Levinson  algorithm, 
LENK  shows  that  the  nonlinear  prediction  errors  can  be 
updated  in  order  through  the  recursion  relation 


mn 

.q 

mn-q 


mn 


0<K    > 

q 


q-i 

mn 

q-1 

mn-q 


(3.20) 


where 


mn 


0<*    >  = 

q      /       mn  2 

(K    ) 


v^ 


mn 


-K 


-  K 


mn 

q 

1 


(3.21) 
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and  the  unique  reflection  coefficients  are  given  by 


mn       q-1   q-1 
K     =  E{e     r    }  (3.22) 

q        mn    mn-q 


[Ref.  12:pp.  146-147] 

A  deficiency  remains  to  be  corrected:  the  goal  is 
to  estimate  y(k),  but  there  is  no  y(k)  term  in  the  two- 
dimensional  data  matrix  of  equation  (3.13).  This  problem  is 
solved  by  adding  another  channel  to  the  lattice  structure 
and  by  exploiting  the  orthogonality  of  the  backward  pre- 
diction errors.  Since  the  backward  prediction  errors  leaving 
the  top  row  of  the  lattice  filter  (Figure  2.9)  are  uncor- 
related,  they  can  be   used  in  a  Fourier  series   to   estimate 

y(k).   These  backward  prediction  errors  can  be  formed  into  a 

2 
length   L  =  (P+l)   vector  defined  as 

X  0     12  L-l     T 

[F        ]  =  [r    ,r     ,r     , . . . , r       ]    (3.23) 
(m,n)-/       00    00-1   00-2       00-L+1 

where  the  subscript  is  in  the  form  mn-q.  Now,  the  error  in 
estimating  y(k)  using  the  data  in  the  Y(k)  tensor  is 
calculated  from 


L  L-l 

e   =  y(k)  -    V     K  y  ,  (3.24) 

k  A'=0   A 


where  the  Fourier  coefficients,  K  /  ,  are  given  by 
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0        -1  L-l 

[K,]  =  [E{y(k)r   },E{y(k)r-    }  ,  .  . . ,E{y(k ) r       }].(3.25) 
A  00  00-1  00-L+1 


The    resulting   generalized   nonlinear    analysis    lattice 

structure  is  depicted  in  Figure  3.3.  [Ref.  12:pp.  147-150] 

The  nonlinear  lattice  structure  will  now  be  examined 

more  closely.    As  can  be  seen  in  Figure  3.3  ,  the  inputs  to 

the   analysis   filter  are  the  normalized  values  of  y(k)   and 

2 
the  (P+l)  +1  ordered  components  of  the  two-dimensional   data 

field  (equation  (3.13)).   The  ordering  of  these  inputs  is  in 

accordance  with  equation  (3.15).    At  a  given  time  k,   these 

2 
(P+l)  +2   inputs   are  evaluated  and  inserted  into   the   left 

side  of  the  lattice  structure.   The  lattice  calculations  are 

conducted  by  starting  at  the   top   row   and   moving  downward 

along  the  northwest-southeast  diagonal,   and  then   advancing 

to   the   top  of  the  next  diagonal  and  so  on.   When   all   the 

lattice   calculations   for  time  k  have  been   completed,   the 

process   is  repeated  for  time  k+1.    Just  as  in  the  case   of 

the  linear  lattice  filter,   the  PARCOR  coefficients  at   each 

lattice   section   act  to  decorrelate  the  two  input   signals. 

The   backward   error  signals  exit  at  the  top  of  the   lattice 

structure,   and  the  forward  error  signals  exit  at  the   right 

side.    The   output  of  interest  for  inverse  filtering  is  the 

forward  error   signal  from  the  top  row  of  the  filter   (i.e., 

the   row  corresponding  to  the  y(k)  input  signal).    This   is 

the  error  arising  from  the  estimation  of  y(k). 


91 


forward  error  signals 


G~0<5 


CD 'CD  CD  Q-G-O  9  0 


■a  -O-O  Q  O  O  O  Q 


tto 


J". 

O 

u 

•a 


0-0-^5-0-0-^ 


o 
£1  -*■ 


0-0-0 


^ 
X 


0-0 


o 


f  c9'9  o  Q 


e 


e-0-0  e  e 


0-0-0-0-0 


I 


c\2 
I 


N 


I 


N 


I 


I 


N 


^ 
X 


W 


N 


X 


N, 


Figure  3.3      Quadratic  Nonlinear  Lattice  Filter 
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3 .   The  Nonlinear  Lattice  Applied  to  Deconvolution 

The  generalized  nonlinear  lattice  filter  can  be 
applied  to  the  inverse  filtering  problem  just  as  the  linear 
lattice  filter  was  in  section  II. D. 6.  The  first  problem  is 
one  of  system  identification  and  parameter  estimation.  Once 
the  model  parameters  have  been  estimated,  they  are  embedded 
in  the  nonlinear  filter  lattice  structure.  Then,  the 
nonlinear  lattice  represents  the  inverse  of  the  system's 
transfer  function.  As  will  be  shown,  the  nonlinear  lattice 
filter  can  model  both  linear  and  nonlinear  systems--  the 
linear  system  is  just  a  special  case  of  the  nonlinear 
system.  This  inverse  filtering  algorithm  was  implemented  by 
the  FORTRAN  programs  NLMAIN,  NLCLAT,  SCHUR,  NORMS,  NLLAT , 
and  URAND .  These  programs  are  listed  in  Appendix  B.  Now, 
this  inverse  filtering  procedure  will  be  described  in  more 
detail . 

In  order  to  identify  the  model  parameters  (i.e.,  the 
PARCOR  coefficients),  the  system  was  excited  with  zero  mean, 
unit  variance  white  noise  sequences.  Both  Guassian  and 
uniform  noise  distributions  were  used  with  good  results. 
One  difficulty  encountered  in  generating  the  nonlinear  data 
was  ensuring  that  the  output  of  the  postulated  nonlinear 
system  remained  bounded  for  the  input  noise  signal.  For  the 
simulations,  the  filter  and  system  were  constrained  to  a 
finite   memory   of   at   most  two   delays   (N=2),   while   the 
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nonlinear  order   was  allowed  to  vary  from  P=0'to  P=4.  (For  N 

greater   than  two,   the   problem  of  ordering  the  elements  of 

the   N-dimensional  Y  tensor   increases  in  complexity.)    The 

system    output     sequence   y(k)   was   then   processed    by 

subroutines    NLCLAT   and   SCHUR,    which    determined    the 

corresponding    autocorrelation    matrix   and   the    partial 

correlation   coefficients,   respectively.    In  an  effort   to 

improve  the  accuracy  of  these  calculations,   noise  sequences 

of  up  to  5,000  points  were  used.   Additionally,   the   PARCOR 

coefficients   were  averaged  over  as  many  as  50   realizations 

of  the  input  noise  random  process.   Through  trial  and  error, 

it   was   found  that  the  best  inverse  filtering  results   were 

obtained   when   the  resulting  reflection   coefficients   were 

truncated  to  two  decimal  places.    The  output  of   subroutine 

2  2 

SCHUR   is  a  ( ( P+ 1 )  +1  x  (P+l)  +1)  upper  triangular  matrix  of 

reflection   coefficients.    This  matrix  is   simply   overlaid 

atop  the  upper  triangular  shaped  nonlinear  lattice  structure 

to   place  the  reflection  coefficients  at  the  correct   filter 

sections . 

With   the   reflection   coefficients   calculated   and 

embedded  in  the  filter,   we  were  able  to  use  the  lattice   in 

an   inverse   filtering   application.   To  recover   the   input 

signal  x(k),   the  system's  output  signal  y(k)  was   processed 

by   subroutines  NORMS  and  NLLAT.    The  function  of  NORMS  was 

2 
to   normalize   the  (P+l)  +1  input   signals  into  the  lattice. 
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Then  subroutine  NLLAT ,  using  the  previously  calculated 
reflection  coefficients,  implemented  the  lattice  structure 
and  carried  out  the  lattice  filter  calculations.  The  forward 
error  signal  out  of  the  top  row  of  the  lattice  yielded  the 
normalized  estimate  of  x(k).  Since  the  inputs  to  the  lat- 
tice filter  are  all  normalized,  the  outputs  must  be 
denormalized .  Therefore,  since  the  input  into  the  top  row  of 
the  lattice  is  divided  by  the  norm  of  y(k),  the  forward 
error  signal  from  the  top  row  of  the  lattice  is  multiplied 
by  the  norm  of  y(k).  (Note  that  if  y(k)  is  a  zero  mean 
sequence,  then  this  norm  is  equivalent  to  its  standard 
deviation.)  Also,  it  was  found  that  to  achieve  a  good 
estimate,  it  was  necessary  to  further  scale  this  error 
signal  by  dividing  it  by  the  norm  of  the  noise  generated 
sequence  y ( k ) . 

In  order  to  verify  the  accuracy  of  the  reflection 
coefficients,  the  noise  generated  output  signal,  y(k),  was 
passed  through  the  inverse  lattice  filter  to  see  how  well 
the  filter  whitened  this  sequence.  The  effectiveness  of  the 
whitening  filter  was  evaluated  by  examining  the  mean-square 
error  ( MSE )  between  the  white  noise  input  signal,  x(k),   and 

A 

the  lattice  filter's  output,  x(k).  The  running  average 
mean-square  error  was  calculated  using  the  equation 


k  a     2 

MSE(k)  ="V/(l/k)  Y.      (x(i)  -  x(i))    .  (3.26) 


i  =  l 
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Since  the  noise  input,  x(k),  has  unit  variance,  the  MSE 
plot  essentially  provides  a  percentage  error  between  x(k) 
and  x(k).  Plots  of  the  MSE  are  included  in  the 
simulation  results.  Having  demonstrated  that  the  lattice 
filter  represented  a  good  inverse  of  the  system  H(z),  the 
system  was  then  driven  by  a  known  input  signal  x(k).  To 
recover  an  approximation  to  this  signal,  the  corresponding 
system  output  was  passed  through  the  lattice  filter,  and  the 
resulting  lattice  filter  output  was  denormalized  and 
rescaled  as  previously  discussed.  Simulation  results  for 
various  systems  are  presented  in  the  following  section. 
4 .   Inverse  Filtering  Simulation  Results 

In  this  section,  the  previous  modeling  and  inverse 
filtering  procedures  are  implemented  and  applied  to  various 
linear  and  nonlinear  systems.  Unless  otherwise  noted,  the 
reflection  coefficients  for  each  of  the  systems  were 
determined  by  using  twenty-five  realizations  (5,000  points 
each  )  of  the  zero  mean,  unit  variance  white  noise  random 
process  as  the  input  excitation  signal.  As  previously 
mentioned,  the  mean-square  error  between  this  white  noise 
input  and  the  lattice  filter's  output  was  plotted  to 
evaluate  the  lattice  filter's  inverse  filtering  performance. 
After  the  system  was  modeled,  it  was  driven  by  the  signal 
x(k)  which  consisted  of  ramps,  pulses,  and  sinusoids  as 
shown  in  Figure  3.4. 
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Figure  3.4    Input  Signal  x(k) 
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a.   System  I:   y(k)  =  x(k)  -  (0.6y(k-l)  +0.08y(k-2)) 

This   is  the  linear  system  first   introduced   in 

Section   II. D. 6.   An  input  Gaussian  white  noise  sequence  was 

used   to  model  the  system.    Figure  3.5   is  a  plot    of   the 

running   average   mean-square  error  between  the  input   noise 

and   output  error  signals.    In  Section  II. D. 6,   the   linear 

lattice  filter's  reflection  coefficients   were   found   to  be 

K   =  -0.555  and   K  =  -0.077.    It  should  be  noted  that  these 

1  2 

values   appear  in  the  first  row  of  the   nonlinear   lattice's 

reflection   coefficient  matrix.   Here,   for  a  first    order, 

nonlinear   lattice  filter,   the  reflection  coefficients   are 

given  by  the  upper  triangular  matrix 


K  = 


0 

0 

-.55 

-.07 

0 

0 

0 

0 

0 

-.43 

0 

0 

0 

-.55 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

The  system  output,  y(k),  corresponding  to  the  input  signal 
of  Figure  3.4  is  shown  in  Figure  3.6.  As  can  be  seen  by  the 
plots  of  x(k)  and  x(k)  in  Figure  3.7,  the  nonlinear  lattice 
filter  did  an  outstanding  job  of  recovering  the  "unknown" 
input  signal  x(k)  from  the  linear  system's  output  y(k). 
b.   System  II:   y(k)  =  x(k)  -  ( 0 . 2y ( k-1 ) y (k-2 ) ) 

This  system  involves  a  "cross-talk" 
nonlinearity ,  that  is,  the  output  is  a  function  of  the 
product   of  two  different  signals.   For   this   system,   both 
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Figure  3.5    System  I  Mean-Square  Error 
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Figure  3.6    System  I  Output  y(k) 
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Figure  3.7    System  I  Comparison  of  x(k)  and  x(k) 
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Gaussian  and  uniform  noise  distributions  were  used  in  the 
modeling  process  in  order  to  compare  the  two  techniques. 
Here,  both  methods  yielded  the  same  reflection  coefficients: 


K  = 


0 

0 

0 

0 

-.21 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

The  only  difference  between  the  two  techniques  was  the  value 
of  the  scaling  factor  (i.e.,  the  norm  of  the  noise  generated 
system  output  y(k)).  The  plot  of  y(k)  is  shown  in  Figure 
3.8.  The  running  average  mean-square  error  and  x(k)  plots 
for  the  Gaussian  noise  derived  model  are  depicted  in  Figures 
3.9  and  3.10,  respectively.  The  corresponding  plots  for  the 
uniform  noise  derived  model  are  given  in  Figures  3.11  and 
3.12.  As  can  be  seen  by  comparing  the  plots,  both  modeling 
variations  lead  to  nearly  identical  results.  Both 
techniques  yielded  excellent  approximations  to  the  input 
signal  x ( k  )  . 

c.   System  III:   y(k)  =  x(k)  -  0 . 2y (k-1 ) y ( k-1 ) 

This  system  involves  a  quadratic  nonlinearity . 
Therefore,  a  second  order  model  was  used.  The  system  output 
would  not  remain  bounded  for  a  Gaussian  noise  input,  so  the 
system  was  modeled  using  a  unit  variance  uniform  noise 
random  process.   The  resulting  reflection  coefficients  are: 
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Figure  3.8    System  II  Output  y(k) 
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Figure  3.9    System  II  Mean-Square  Error 
(Gaussian  Noise  Model) 
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Figure  3.10   System  II  Comparison  of  x(k)  and  x(k) 
(Gaussian  Noise  Model) 
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Figure  3.11    System  II  Mean-Square  Error 
(Uniform  Noise  Model) 
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Figure  3.12    System  II  Comparison  of  x(k)  and  x(k 

(Uniform  Noise  Model) 
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K  = 


0 

-.20 

.08 

0 

-.04 

-.18 

0 

0 

0 

0 

0 

0 

-.20 

-.18 

0 

.69 

.54 

.  17 

.34 

-.42 

0 

0 

0 

.12 

-.22 

-.39 

-.04 

.01 

.64 

.31 

0 

0 

0 

0 

-.38 

-.09 

-.37 

.60 

-.15 

.32 

0 

0 

0 

0 

0 

.27 

.15 

-.67 

-.25 

-.  13 

0 

0 

0 

0 

0 

0 

.57 

-.  17 

-.37 

.49 

0 

0 

0 

0 

0 

0 

0 

-.36 

-.37 

.50 

0 

0 

0 

0 

0 

0 

0 

0 

.51 
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Figure  3.13  provides  the  running  average  mean-square  error 
plot.  The  output  y(k)  of  the  nonlinear  system  is  shown  in 
Figure  3.14.  Figure  3.15  depicts  the  comparison  between  the 
system  input  x(k)  and  the  lattice  filter  output  x(k).  The 
x(k)  curve  has  the  same  shape  as  x(k),  however,  it  is 
slightly  offset.  In  an  attempt  to  improve  the  approximation 
process,  the  number  of  realizations  of  the  noise  random 
process  used  to  model  the  system  was  doubled  from  twenty- 
five  to  fifty  and  the  simulation  was  repeated.  Although 
several  reflection  coefficient  values  changed  by  .01,  there 
was  no  perceptible  improvement  in  the  estimation  of  x(k). 

d.   System  IV:  y(k)=x(k)  -  ( 0 . 5+0 . 6y ( k-1 ) + . 08y ( k-2 ) ) 

This   example   consists  of  the  linear   System   I 

modified  by  the  addidtion  of  a  constant  bias  term.   Gaussian 

noise   was   used   to   model   the   system.     The   reflection 

coefficients  determined  for  the  first  order  model  are: 
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Figure    3.13         System    III    Mean-Square    Error 


109 


it 


Figure  3.14    System  III  Output  y(k) 
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Figure  3.15    System  III  Comparison  of  x(k)  and  x(k) 
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K  = 


0 

-.23 

-.55 

-.07 

0 

0 

0 

-.23 

-.40 

-.43 

0 

0 

0 

-.46 

.01 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

This  is  basically  the  same  matrix  as  obtained  for  System  I, 
but  modified  by  the  addition  of  several  terms  which  attempt 
to  model  the  constant  offset.  (Note  the  reappearance  of  the 
values  of  -.55  and  -.07  in  the  top  row  of  the  lattice 
matrix.)  The  running  average  mean-square  error  is  shown  in 
Figure  3.16.  The  system  output  y(k)  is  plotted  in  Figure 
3.17.  Figure  3.18  is  a  plot  of  the  lattice  output  x(k) 
and  the  system  input  x(k).  The  x(k)  curve  is  an  excellent 
replica  of  the  original  input  signal,  but  it  is  offset  by  a 
constant  value.  Although  the  small  mean-square  error  evident 
in  Figure  3.16  indicates  that  the  lattice  is  a  good  inverse 
filter,  the  lattice  was  unable  to  completely  remove  the  bias 
term.  Increasing  the  order  of  the  model  (i.e.,  overmodeling 
the  system)  did  not   improve  the  results. 

e.  System  V:  y(k)=x(k)  -  ( 5 . 0+0 . 6y ( k-1 ) + . 08y ( k-2 ) ) 
In  order  to  further  investigate  the  constant 
offset  nonlinearity ,  the  example  of  System  IV  was  repeated 
using  a  bias  of  5.0  instead  of  0.5.  Again,  Guassian  white 
noise  was  used  in  determining  the  model's  parameters.  The 
reflection  coefficients  of  the  first  order  nonlinear  lattice 
are  given  by: 
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Figure  3.16    System  IV  Mean-Square  Error 
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Figure  3.17    System  IV  Output  y(k) 
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Again,  the  reflection  coefficient  values  of  -.55  and  -.07 
occur  in  the  top  row  of  the  matrix  as  these  elements 
apparently  model  the  linear  portion  of  the  sytem.  The 
running  average  mean-square  error  is  displayed  in  Figure 
3.19,  and  the  output  of  the  nonlinear  system  is  shown  in 
Figure  3.20.  Figure  3.21  provides  the  comparison  between  the 
system  input  and  the  output  of  the  inverse  filter.  As  with 
System  IV,  the  x(k)  curve  produced  by  the  deconvolution 
process  has  the  same  shape  as  x(k),  but  is  offset  from  the 
desired  curve  by  a  small  constant  value.  Here,  the  system 
bias  is  5.0,  and  the  x(k)  and  x(k)  curves  differ  by  about 
0.5.  This  is  a  relative  improvement  over  the  results 
obtained  for  System  IV.  System  IV  had  a  constant  bias  of 
only  0.5  but  the  x(k)  curve  was  offset  from  the  x(k)  curve 
by  approximately  0.4. 
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Figure    3.19         System   V   Mean-Square    Error 
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Figure  3.20    System  V  Output  y(k) 
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Figure  3.21    System  V  Comparison  of  x(k)  and  x(k) 
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IV.   CONCLUSION 

In  this  thesis,  several  linear  least-squares 
deconvolution ,  or  inverse  filtering,  techniques  were 
reviewed.  Particular  emphasis  was  placed  on  the  lattice 
filter.  It  was  shown  that  when  a  system  is  defined  by  an 
autoregressive ,  linear  time-invariant  model,  this  model  can 
be  transformed  into  an  equivalent  lattice  filter 
representation.  Furthermore,  the  resulting  analysis  lattice 
filter  acts  as  a  whitening  filter,  and  is  the  inverse  of  the 
system's  autoregressive  transfer  function.  An  example 
demonstrated  the  effectiveness  of  the  lattice  filter  in  a 
linear  deconvolution  application.  Unfortunately,  the  linear 
lattice  filter  is  unable  to  model  nonlinear  systems  and, 
therefore,  is  not  a  viable  nonlinear  inverse  filtering 
technique . 

The  discussion  of  the  linear  lattice  filter  led  to  the 
development  of  the  generalized  lattice  filter,  which  in  turn 
led  to  the  derivation  of  the  nonlinear  lattice  filter. 
The  goal  of  this  thesis  was  to  implement  the  nonlinear 
lattice  filter,  and  then  apply  it  to  the  linear  and 
nonlinear  deconvolution  problem.  This  was  accomplished  with 
generally  very  good  results.  It  was  shown  that  the 
nonlinear    lattice    filter   was   suitable    for    modeling 
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discrete   autoregressive   systems   where   the   output    y(n) 

consisted   of   a   weighted   sum   of   the   cross-products   of 

P  q 

the  terms  y  (n-1)  and  y  (n-2),   where  the  powers  p  and  q  can 

take   on  the  values   {0,1,2,3,4}.    Examples  were   presented 

demonstrating  the  ability  of  the  nonlinear  lattice  filter  to 

effectively   act   as  an  inverse  filter  for  both   linear   and 

nonlinear  autoregressive  systems. 
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APPENDIX  A 
LINEAR  LATTICE  FILTER  FORTRAN  PROGRAMS 


p  ****************************************************************** 

C  *  * 

C  *   LT  SCOT  L  JOHNSON  * 

C  *   LINEAR  INVERSE  * 

C  *   06  APRIL  1986 

C 

C 


*  * 

****************************************************************** 

THE  PURPOSE  OF  THIS  PROGRAM  IS  TO  MODEL  A  LINEAR  AUTOREGRESSIVE 

C  SYSTEM  BY  THE  EQUIVALENT  ANALYSIS  LATTICE  FILTER. 

C  THE  SYSTEM'S  TRANSFER  FUNCTION  IS  GIVEN  BY  H(Z)  =  1/A(Z) 

C  WHERE  A(Z)  =  (Z**2  +  A1*Z  +  A2)/Z**2. 

C         THE  LATTICE  FILTER  REFLECTION  COEFFICIENTS  ARE  DETERMINED 

C  BY  DRIVING  H(Z)  WITH  ZERO  MEAN,  UNITY  VARIANCE  GAUSSIAN  WHITE 

C  NOISE.  THE  OUTPUT  IS  PROCESSED  BY  SUBROUTINE  ATOCOR  WHICH  COMPUTES 

C  THE  COMPONENTS  OF  THE  SAMPLED  AUTOCORRELATION  MATRIX  R.  SUBROUTINE 

C  LEVIN  IMPLEMENTS  THE  LEVINSON  ALGORITHM  AND  EVALUATES  THE  LATTICE 

C  FILTER  COEFFICIENTS  FROM  THE  AUTOCORRELATION  MATRIX.  FINALLY, 

C  SUBROUTINE  LATICE  IMPLEMENTS  THE  ANALYSIS  STRUCTURE  OF  THE 

C  LATTICE  FILTER  WHICH  IS  EQUIVALENT  TO  THE  INVERSE  FILTER  H(Z). 


C  NOW  WHEN  H(Z)  IS  DRIVEN  BY  AN  UNKNOWN  SIGNAL  XfN],  THE 

C  SYSTEM  OUTPUT  Y(N)  CAN  BE  FED  INTO  THE  PREVIOUSLY  DETERMINED 

C  LATTICE  FILTER;  SINCE  THE  LATTICE  FILTER  IS  THE  INVERSE  OF  H(Z) 

C  THE  LATTICE  FILTER  OUTPUT  SHOULD  BE  A  GOOD  ESTIMATE  OF  X(N). 

C 

C  ***VARIABLE  DEFINITIONS**** 

C  Al.  A2  =  COEFFICIENTS  OF  THE  PREDICTION  ERROR  POLYNOMIAL  A(Z] 

C  DSEED,DSEED1  =  SEED  VALUES  USED  BY  THE  IMSL  WHITE  GAUSSIAN  NOISE 

C  GENERATOR  FUNCTION  GGNQG 

C  E  =  VECTOR  OF  MEAN-SQUARED  PREDICTION  ERRORS  E(0) ,E(1) ,. . . E(ORDER) 

C  ESUM  =  VECTOR  USED  IN  DETERMINING  AVERAGE  E 

C  GAMMA  =  VECTOR  OF  LATTICE  REFLECTION  COEFFICIENTS.  GAMMA(I)  IS 

C  THE  REFLECTION  COEFFICIENT  OF  THE  I'TH  LATTICE  STAGE. 

C  GAMSUM  =  VECTOR  USED  IN  DETERMING  AVERAGE  GAMMA 

C  GRAFP  =  REAL  VALUE  WHICH  DEFINES  LENGTH  OF  X-AXIS  FOR  PLOTTING 

C  L  =  LOWER-TRIANGULAR  MATRIX  WHOSE  ROWS  ARE  THE  REVERSE  OF  ALL  THE 

C  PREDICTION  ERROR  FILTERS  FROM  ORDER  ZERO  TO  THE  HIGHEST  ORDER 

C  LSUM  =  MATRIX  USED  IN  DETERMINING  THE  AVERAGE  L  MATRIX 

C  MSE  =  MEAN-SQUARE  ERROR  BETWEEN  X(N)  AND  XHAT(N) 

C  MSEMAX  =  MAXIMUM  VALUE  OF  MSE 

C  MSESUM  =  RUNNING  SUM  USED  IN  CALCULATING  MSE 

C  NINDEX  =  ARRAY  OF  REAL  NUMBERS  USED  IN  DISSPLA  PLOTTING  ROUTINES 

C  NOISE  =  ARRAY  OF  MEASUREMENT  NOISE  ADDED  TO  OUTPUT  OF  H(Z) 

C  NUMPTS  =  NUMBER  OF  POINTS  IN  INPUT  NOISE  SEQUENCE 

C  ORDER  =  ORDER  OF  THE  LATTICE  FILTER 

C  ORDERP  =  ORDER  +  1 

C  PLTPTS  =  NUMBER  OF  POINTS  USED  IN  PLOTTING  ROUTINES 

C  R  =  VECTOR  OF  AUTOCORRELATION  LAGS  R(0),R(1} .. . . R(ORDER) 

C  RMAX  =  MAXIMUM  MAGNITUDE  OF  ELEMENTS  OF  Rl  l6  BE  USED  IN  PLOTTING 

C  STDDEV  =  STANDARD  DEVIATION  OF  THE  MEASUREMENT  NOISE 

C  TRIAL  =  NUMBER  OF  REALIZATIONS  OF  THE  WHITE  GAUSSIAN  NOISE  RANDOM 

C  PROCESS  USED  IN  MODELING  THE  SYSTEM  H(Z) 

C  X  =   INPUT  SEQUENCE  INTO  THE  SYSTEM  H(Z) 

C  XMAX,  XMIN  =  RANGE  OF  X  VALUES  USED  IN  DISSPLA  PLOTTING  ROUTINES 

C  XHAT  =  ESTIMATE  OF  X-  OUTPUT  OF  THE  ANALYSIS  LATTICE  FILTER  A(Z) 

C  Y  =  OUTPUT  SEQUENCE  OF  H(Z) 

C  YMAX,  YMIN  =  RANGE  OF  Y  VALUES  USED  IN  DISSPLA  PLOTTING  ROUTINES 

C 

C  ***VARIABLE  DECLARATIONS**** 

INTEGER  I, N,K, NUMPTS, ORDER, ORDERP, PLTPTS, TRIAL 
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X( 5000). Y( 5000), Al.A2.XMIN,XMAX,YMAX.YMIN,NINDEXf5( 
GAMMA( i6),GAMS0M( l6) . ft( 10)  L( 10  10) AE( 10) ,RMAX ,GRAf 
XHATL5000)fLSUMClO,ld)lESOfl(10)  NOiSE(5006;,ST6DEV 

MSE(5000)>1SESUM,M$EPLt 


REALM  X( 5000). Y( 5000), Al.A2.XMIN,XMAX,YMAX.YMIN,NINDEXf 5000) 
REALM  GAMMA(i6),GAMSU^l6),^10)!LC10!l0),E(10)!RMAXJGBAFP 

REALM 

REAL*8  DSEED^DSEEbl' 
C 
C 
C     DEFINE  THE  LENGTH  OF  THE  TIME  SEQUENCES 

NUMPTS  =  5000 

PLTPTS  =  2000 

GRAFP  =  FLOAT(PLTPTS+l) 

C     DEFINE  THE  ORDER  OF  THE  LATTICE  FILTER  AND  THE  A(Z)  COEFFICIENTS 

ORDER=2 

ORDERPORDER+1 

Al  =  0.6 

A2  =  0.08 
C     INITIALIZE  VARIABLES 

DO  6  I=l,ORDERP 


GAMMAU)  =  0. 
GAMSUM(l)  =  ( 
ESUM(I)  =  0. 


RCI)  =  0. 

DO  4  K=l,ORDERP 
LSUM(I.K)  =  0. 
4         CONTINUE 
6     CONTINUE 
C 

XMAX  =  0. 
YMAX  =  0. 
MSESUM  =  0. 
MSEMAX  =  0. 
C     DEFINE  THE  SEEDS  FOR  THE  NOISE  SIGNALS,  THE  NUMBER  OF  NOISE 
C     REALIZATIONS  USED  IN  MODELING  H(Z),  AND  THE  STANDARD  DEVIATION  OF 
C     THE  MEASUREMENT  NOISE 
DSEED  =  1243073. 5D0 
DSEED1  =  724389. 4D0 
TRIAL  =  25 
STDDEV  =  0.05 
C 

DO  25  I=1,TRIAL 

DSEEd  =  DSEED/DFLOAT(I) 
DSEED1  =  DSEED1/DFL0AT(I) 
DO  10  K=l, NUMPTS 

X(K)  =  GGNQF(DSEED) 

10  CONTINUE 

DO  11  K=l. NUMPTS 

NOIS£(K)  =  STDDEV  *  GGNQF(DSEEDl) 

11  CONTINUE 

Y(l)  =  X(l)  +  NOISE(l)  ,  x 
Y(2)  =  X  2)  -  A1*Y(1)  +  NOISE(2) 
DO  12  K  =  3, NUMPTS 

Y(K)  =  X(K)  -  (A1*Y(K-1)  +  A2*Y(K-2))  +  NOISE(K) 

12  CONTINUE 

C        DETERMINE  THE  AUTOCORRELATION  VECTOR  OF  THE  SEQUENCE  Y(N) 

CALL  ATOCOR(NUMPTS,Y, ORDER, R.RMAX) 
C        DETERMINE  THE  L  MATRIX  AND  E  VECTOR 
CALL  LEVIN(ORDER,R,L,E) 
DO  18  K=l. ORDER 

GAMMACK)=-1*L(K+1.1) 
GAMSUM(K)  =  GAMSUM(K)  +  GAMMA(K) 
18       CONTINUE 
C 

DO  22  N=1.0RDERP 

ESUMfN)  =  ESUM(N)  +  E(N) 
DO  20  K=l,ORDERP  /M       lllll(, 

LSUM(N,K)  =  LSUM(N,K)  +  L(N,K) 
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20  CONTINUE 

22      CONTINUE 

25  CONTINUE 
C 

C     CALCULATE  THE  AVERAGE  VALUES  OF  E,  L,  AND  GAMMA 
DO  28  N=l,ORDERP 

E(N)  =  ESUM(N)/FLOAT(TRIAL) 
DO  26  K=l,ORDERP 

L(N,K)  =  LSUM(N,K)/FLOAT(TRIAL) 

26  CONTINUE 
28    CONTINUE 

DO  30  K=l. ORDER 

GAMMX(K)  =  GAMSUM(K)/FLOAT(TRIAL) 
30    CONTINUE 
C 

C    WITH  INPUT  Y  AND  LATTICE  PARAMETERS  GAMMA4  DETERMINE  THE 
C       OUTPUT  OF  THE  ANALYSIS  MODEL  LATTICE  FILTER 
C     CALL  LATICE(NUMPTS, ORDER, GAMMA, Y,XHAT) 

C         PRINT  RESULTS 

C 

C     PRINT  THE  R.E,  AND  GAMMA  VECTORS 

WRITE(8,35) 
35  FORMATm,1!1  ,T4,'N"  ,T10,'R(N)'  ,T20,'E(N)'  ,T30,  'GAMMA(N)  ■ ) 

NULL=D 

WRITE(8,40)  NULL.R(l).EjflJ 
40    FORMATCTl . ' 0 ' . T4 . f 2 ,Tld ,  3F10.  4) 

DO  45  K=2,ORDERP 
N=K-1 

WRITE(8,40)  N,R(K),E(K),GAMMA(N) 
45    CONTINUE 
C 
C     PRINT  THE  L(I,J)  MATRIX 

WRITE(8V48) 

48         FORMAT  CTl/l'.'LCl,  J  )=') 

DO  55  I=l,ORDERP 

WRltEi8,50)  (L(I,J),J=1.0RDERP) 
50  F0RMAT(tl,r0^I0CFi0.4,4XJ) 

55    CONTINUE 

C 

C 

C     DO  THE  DECONVOLUTION  PROBLEM.  FIRST  DEFINE  THE  DETERMINISTIC 

C     INPUT  X(N),  THEN  DETERMINE  THE  SYSTEM  OUTPUT. 

DO  85  K=1,PLTPTS 

X(K)  =  1.0  *  SIN(0.  0126*FLOAT(K)) 
85    CONTINUE 
C 

Y(l)  =  X(l)  +  NOISE(l) 

Yf2]  =  Xf2)  -  A1*Y(1)  +  NOISE(2) 

DO  $0  K=3,PLTPTS 

Y(K)  =  X(K)  -  (A1*Y(K-1)  +  A2*Y(K-2))  +  NOISE(K) 
90    CONTINUE 
C 
C     INPUT  Y(N)  INTO  THE  LATTICE  TO  RECOVER  THE  ESTIMATE  OF  X(N) 

CALL  LAT1CE(PLTPTS, ORDER, GAMMA, Y,XHAT) 

C     CALCULATE  THE  MEAN-SQUARE  ERROR  (MSE)  BETWEEN  X(N)  AND  XHAT(N) 
DO  93  K=1,PLTPTS 

MSE$UM  =  fX[K)-XHAT(K))*(X(K)-XHAT(K))  +  MSESUM 
MSE(K)  =  SQRTfMSESUM/FLOATCK)) 
IF(MSE(K).GT.  MSEMAX)  MSEMAX  =  MSE(K) 

93    CONTINUE 

C 

C     DEFINE  PLOTTING  PARAMETERS 
DO  95  K=  l.PLTPTS 

NINDEX(K)  =  FLOAT(K-l) 
IF(ABS(X(K)).GT.  XMAX)  XMAX  =  ABS(X(K)) 
IF(ABS(XHAT(K)).GT. XMAX)  XMAX  =  ABS(XHAT(K)) 
IF(ABS(Y(K)).GT.  YMAX)  YMAX  =  ABS(Y(K)) 
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95    CONTINUE 

XMIN  =  -1.0  *  XMAX 

YMIN  =  -1.0  *  YMAX 
C 

C     ****  DISSPLA  PLOTTING  ROUTINES  **** 
C 

C     PLOT  THE  SYSTEM  INPUT  AND  THE  LATTICE  OUTPUT 
C     CALL  TEK618 

CALL  COMPRS 

CALL  RESETpALL') 

CALL  PAGE(8.  0,6.  5) 

CALL  XINTAX 

CALL  AREA2D(6.  75,5.0) 

CALL  XNAMEpTIME  SAMPLES$' ,100) 

CALL  YNAME('X(N)$r,100) 

CALL  CROSS 

CALL  GRAFf 0. ,500. ,GRAFP, XMIN, 'SCALE' ,XMAX) 

CALL  CURVEfNlNDEX  X.PLTPTS.O) 

CALL  ENDPL(O)  . 

C 
C     PLOT  THE  SYSTEM  OUTPUT  Y(N) 

CALL  NOBRDR 

CALL  AREA2D(6.  75,5.0) 

CALL  XNAMEpTIME  SAMPLES$' ,100) 

CALL  YNAMEf 'Y(N)$',100) 


2) 


C     CALL  HEADINfYCN)  =  X(N)  -  (0.  6*Y(N-l)+0.  08*YCN-2)1$'  .100.1.  5.11 

CALL  HEADIN('Y(N)=X(N)-[0.6*YrN-l)+0.08',(Y(N-2))+N6lSE$,,!L00,1.5, 

CALL  HEADIN( ' NOISE  =  N(5,0.  0025)$\  100, 1.5, 2) 

CALI  CROSS 

CALL  GRAFfo.  ,500. ,GRAFPAYMIN, 'SCALE' ,YMAX) 

CALL  CURVEC NlNDEX I Y , PLT^TS , 0) 

CALL  ENDPL(O) 
C 
C     PLOT  LATTICE  FILTER  OUTPUT,  XHAT(N) 

CALL  AREA2D(6.  75,5.0) 

CALL  XNAMEpTIME  SAMPLES* '  ,100) 

CALL  YNAME('XHAT(N)$', 100) 

CALL  CROSS 

CALL  GRAF(0.  ,500.  .GRAFP. XMIN. "SCALE1 ,XMAX) 

CALL  CURVEfNlNDEX  XHAT,PLTPT$,0) 

CALL  ENDPL(O) 
C 
C     PLOT  THE  MEAN-SQUARE  ERROR  BETWEEN  X(N)  AND  XHAT(N) 

CALL  AREA2D(6.  75,5.0) 

CALL  XNAMEPTIME  SAMPLES?',  100) 

CALL  YNAMEp MSE( N)$',  100) 

CALI    CROSS 
C  CALL  GRAF(0. ,500. , GRAFP .0. ,' SCALE' .MSEMAX) 

CALL  GRAFfO!  ,500!  ,GRAFP,o!  .  ■  SCALE^  ,0.  08) 

CALL  CURVE(NlNDEX,MSE,PLTPtS,0) 

CALL  ENDPL(O) 

CALL  DONEPL 
C 

STOP 

END 
C 
C 

C 

C  SUBROUTINE  ATOCOR 

C 

C    GIVEN  A  TIME  SEQUENCE  Y(N),  THIS  PROGRAM  CALCULATES  THE  SAMPLED 

C    AUTOCORRELATION  MATRIX  TERMS  R(0) ,R(1) ,. . .  ,R(ORDER). 

C    WRITTEN  06  APRIL  1986 
C 

p*********************  **********************************************  A*** 

SUBROUTINE  ATOCOR(NUMPTS,Y, ORDER, R.RMAX) 
C     VARIABLE  DEFINITIONS: 
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C  Y(5000)=  INPUT  SEQUENCE 

C  Rf  10)  =  AUTOCORRELATION  MATRIX 

C  NUMPTS  =  NUMBER  OF  POINTS  IN  THE  SEQUENCE  Y(N) 

C  ORDER  =  MAXIMUM  LAG  FOR  WHICH  R  IS  EVALUATED 

C  ORDERP  =  0RDER+1:THE  LENGTH  OF  THE  R  VECTOR 

C  SUM=  RUNNING  TOTAL  OF  PRODUCTS  FOR  A  GIVEN  LAG 

C  RMAX=  MAXIMUM  VALUE  OF  R;  USED  FOR  DISSPLA  PLOTING 


C 


INTEGER  NUMPTS, I ,J,K,N, ORDER, ORDERP 

\R(i6X 


REAL  Y(5000),R(lO),SUM,RMAX 
RMAX=0. 
0RDERP=0RDER+1 
DO  75  K=l, ORDERP 
SUM=0. 

N=NUMPTS-K+1 
DO  70  J=1,N 

SUM=SUM+(Y(K+J-1)*Y(J)) 
70        CONTINUE 

R(K)=SUM/NUMPTS 

IF(ABS(R(K)).  GT.  RMAX)  RMAX=ABS(R(K)) 
75    CONTINUE 
C 

RETURN 
END 
C 

C 

p*****  ***********************************************  ******************* 

C 

C     SUBROUTINE  LEVIN 

C 

C     THIS  SUBROUTINE  IMPLEMENTS  LEVINSON'S  ALGORITHM.  IT  GENERATES  ALL 

C     THE  PREDICTION  ERROR  FILTERS  UP  TO  A  GIVEN  ORDER,  FROM  THE 

C     AUTOCORRELATION  LAGS. 

C 

C     BASED  ON  A  PROGRAM  WRITTEN  BY  S.J.  ORFANIDIS  (REF.  1:  P.  333) 

C 

C    MODIFIED  06  APRIL  1985 

C 

p************* ******************************************* *************** 

SUBROUTINE  LEVIN(ORDER,R,L,E) 
C     ***VARIABLE  DEFINITIONS'**'** 
C 

C    ORDER  =  ORDER  OF  LATTICE  FILTER 

C     R  =  VECTOR  OF  AUTOCORRELATION  LAGS  Rf 0) Ml) .. . . ,R(ORDER) 
C     L  =  UNIT  LOWER-TRIANGULAR  MATRIX.   ITS  i-TH  ftOW  HOLDS  THE  I-TH 
C         PREDICTION-ERROR  FILTER  IN  REVERSE  ORDER.  ITS  FIRST  COLUMN 
C         HOLDS  THE  NEGATIVES  OF  THE  REFLECTION  COEFFICIENTS, 
C         GAMMA(I)  =  -L(I+1,1)  FOR  1=1,2.. .. .ORDER 
C     E  =  VECTOR  OF  PREDICTION  ERRORS  E( 0) ,E(i],. .  .  ,E(ORDER) 
C     THE  MATRIX  L  AND  THE  DIAGONAL  MATRIX  D  FORMED  BY  THE  E'S  DEFINE 
C     A  UL  CHOLESKY  FACTORIZATION  OF  THE  INVERSE  OF  THE  AUTOCORRELATION 
C     MATRIX:   R  INVERSE  =  L  TRANSPOSE  *  D  INVERSE  *  L 
C 
C     ***VARIABLE  DECLARATIONS**** 

REAL  R(10),E(10),L(10A10)AGAP,GAMMA 

INTEGER  rlPLUS.lMlNUS, J, dRDEA, ORDERP 

orderp=or6er+i 
c   set  the  upper  triangle  of  the  l  matrix  to  zero 

DO  60  1=1, ORDER 
IPLUS=I+1 
IMINUS=I-1 

DO  60  J=IPLUS.ORDERP 
L(I,J)=0\ 
60    CONTINUE 


C 


L(l,l)=l. 

L(2  2)=1.0 

L(2  1]=-R(2)/R(1) 

E Cl)=R(l) 
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GAP=0. 

IMINUS=I-1 

DO  63  K=1,IMINUS 

GAP=GAP+R(K+1)*L(I-1,K) 

63  CONTINUE 

GAMMA=GAP/E(I-1) 
L(I,1)=-1*GAMMA 
DO  64  K=2,IMINUS 

L(I  K)=L(I-1,K-1)-GAMMA*L(I-1,I-K) 

64  CONTINUE 


L(1. 11=1.0 


■1)*(1-GAMMA**2) 
68    CONTINUE 


RETURN 
END 


r>  *********************************************************** ************ 

c 

C  SUBROUTINE  LATICE 

C 

C  THIS  PROGRAM  IMPLEMENTS  A  SINGLE  CHANNEL  LATTICE  STRUCTURE. 

C  WHEN  GIVEN  THE'  LATTICE  COEFFICIENTS  AND  THE  INPUT  SEQUENCE, 

C  THE  PROGRAM  DETERMINES  THE  OUTPUT  SEQUENCE. 

C    WRITTEN  06  APRIL  1986 

C 

p*********************************************************************** 

c 

SUBROUTINE  LATICEfNUMPTS, ORDER, GAMMA, Y, OUTPUT) 

C  ***VARIABLE  DEFINITIONS**5** 

C  NUMPTS=  NUMBER  OF  POINTS  IN  THE  SEQUENCES;  MAX  IS  5000 

C  ORDER=  ORDER  OF  THE  LATTICE-  MAX  IS  9 

C  GAMMAfORDER)=  LATTICE  COEFFlCENT  ARRAY 

C  F  =  FORWARD  ERROR 

C  B  =  BACKWARD  ERROR 

C  DELAY(ORDER)=  ARRAY  OF  DELAYED  BACKWARD  ERROR  SIGNALS 

C  TEMP(ORDER)  =  ARRAY  WHICH  TEMPORARILY  HOLDS  THE  BACKWARD  ERROR 

C  Y(NUMPTS)=INPUT  DATA  ARRAY 

C  O0TPUT(NUMPTS)=ARRAY  OF  LATTICE  OUTPUT  DATA 

C     ***VARIABLE  DECLARATIONS**** 

INTEGER  NUMPTS, ORDER, I.K.M         ■   ,  , 

REAL  GAMMAflO)  F,B,DELAY(10),TEMP(10),Y(5000),OUTPUT(5000) 
C     INITIALIZE  ARRAYS 

DO  80  1=1 .ORDER 
DELAy(I)=0. 
TEMP(I)=0. 
80    CONTINUE 
C 
C     DO  TIME  ITERATION 

DO  88  K=l, NUMPTS 

f=y£io 

c      for~each  time  instant, recursively  increase  the  lattice  order 

DO  85  M=l, ORDER 
TEMlVM)=B 


B=DELAY(M)-(GAMMA(M)*F) 
F=F:tGAMMA(M)*DELAY(M) ) 


DELAY(M)=TEMP(M) 
85        CONTINUE 

OUTPUT(K)=F 
88    CONTINUE 
RETURN 
END 
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c 


APPENDIX  B 
NONLINEAR  LATTICE  FILTER  FORTRAN  PROGRAMS 

********************************************************************* 

* 

C  PROGRAM  NLMAIN  * 

C 

C  THIS  PROGRAM  DEFINES  THE  SYSTEM  PARAMETERS  AND  THE  SYSTEM'S  * 

C  INPUTS  AND  OUTPUTS.   IT  CALLS  SUBROUTINES  TO  DETERMINE  THE      * 

C  THE  CORRESPONDING  NONLINEAR  ANALYSIS  LATTICE  MODEL.            * 

C  * 

C  WRITTEN  30  APRIL  1986 

C  * 

p************ ***************************************************** ****** 

C 

C  THIS  PROGRAM  INPUTS  A  WHITE  NOISE  SEQUENCE  X(K)  INTO  A  AUTOREGRES- 

C  SIVE^  NONLINEAR  SYSTEM  H(Z)  .  THE  SYSTEMYS  OUTPUT,  Y(K),IS 

C  PROCESSED  TO  DETERMINE  THE  AUTOCORRELATION  MATRIX  AND  THE  NONLINEA 

C  LATTICE'S  REFLECTION  COEFFICIENT  MATRIX.  SINCE  THE  OUTPUT  OF  THE 

C  LATTICE  IS  WHITE  NOISE,  THE  LATTICE  IS  EQUIVALENT  TO  THE  INVERSE 

C  OF  THE  SYSTEM. 

C 

C  ***  VARIABLE  DECLARATIONS**** 

C  DSEED  =  SEED  USED  BY  THE  IMSL  WHITE  GAUSSIAN  NOISE  FUNCTION  GGNQF 

C  GRAFP  =  MAX  X-AXIS  VALUE  FOR  X  AND  Y  GRAPHS 

C  GRAFN  =  MAX  X-AXIS  VALUE  FOR  MSE  GRAPH 

C  H  =  ARRAY  OF  AUTOREGRESSIVE  PARAMETERS  THAT  DEFINE  THE  SYSTEM 

C  IY  =  SEED  USED  BY  THE  UNIFORM  WHITE  NOISE  FUNCTION  URAND 

C  MN  =  N  *  N 

C  MNP1  =  MN  +  1;  NUMBER  OF  ROWS  IN  THE  NONLINEAR  LATTICE 

C  MSE  =  MEAN  SQUARE  ERROR  BETWEEN  THE  INPUT  NOISE  SIGNAL  AND  THE 

C  FORWARD  ERROR  SIGNAL  FROM  THE  TOP  ROW  OF  THE  LATTICE 

C  MSESUM  =  RUNNING  TOTAL  USED  IN  CALCULATION  MSE 

C  MSEMAX  =  MAXIMAUM  VALUE  OF  MSE  FOR  USE  IN  PLOTTING 

C  MSEPLT  =  REAL*4  VALUES  OF  MSE  USED  IN  DISSPLA  PLOTTING 

C  NORM  =  ARRAY  OF  FACTORS  THAT  NORMALIZE  THE  LATTICE  INPUTS 

C  N  =  DIMENSION  OF  THE  SQUARE  Y  TENSOR  MATRIX 

C  NUMPTS  =  NUMBER  OF  POINTS  USED  IN  CALCULATING  THE  RHO  MATRICES 

C  NINDEX  =  TIME  INDEX  ARRAY  USED  IN  DISSPLA  PLOTS 

C  PLTPTS  =  NUMBER  OF  DATA  POINTS  USED  IN  DISSPLA  PLOTS 

C  RANGE  =  +/-  RANGE  OF  UNIFORM  WHITE  NOISE 

C  RHO  =  ARRAY  OF  REFLECTION  COEFFICIENTS 

C  RHOSUM  =  ARRAY  USED  IN  CALCULATING  THE  AVERAGE  RHO  MATRIX 

C  SCALE  =  RECIPROCAL  OF  THE  NORM  OF  THE  SYSTEM  OUTPUT  FOR  WHEN 

C  THE  SYSTEM  IS  EXCITED  BY  WHITE  NOISE 

C  STDDEV  =  STANDARD  DEVIATION  OF  GUASSIAN  WHITE  NOISE 

C  TRIAL  =  NUMBER  OF  NOISE  REALIZATIONS  USED  IN  CALCULATING  AVG  RHO 

C  X  =  INPUT  INTO  SYSTEM  H(Z) 

C  XHAT  =  LATTICE  OUTPUT  THAT  APPROXIMATES  X 

C  XHPLOT  =  ARRAY  OF  REAL*4  VALUES  OF  XHAT  USED  IN  DISSPLA  PLOTS 

C  XMAX.XMIN  =  RANGE  OF  X  AND  XHPLOT  VALUES  USED  IN  DISSPLA  PLOTS 

C  Y  =  6UTPUT  OF  HCZ) 

C  YPLOT  =  ARRAY  OF  REAL*4  VALUES  OF  Y  USED  IN  DISSPLA  PLOTS 

C  YMAX,YMIN  =  RANGE  OF  YPLOT  VALUES  USED  IN  DISSPLA  PLOTS 

C  YSUM  =  RUNNING  TOTAL  USED  IN  CALCULATING  THE  SYSTEM  OUTPUT 

C 

C  ***VARIABLE  DECLARATIONS**** 

INTEGER  TRIAL, NUMPTS, PLTPTS, IY.IM1.JM1 

REAL*4  X(5000 j, MSEMAX, RANGE! STDDEV 'YMAX.YMIN.XMAX.XMIN.GRAFN 
REAL*4  XHPLOT(^OOO)  ^iNDEX^OOOhM^EPLTtBOOOl^YPLdTrBOOO) , GRAFP 
PEAL*8  Y(5000)AR( 26! 26)  ARHQSUM(26, 26),  XHAT(  5000], NORM(26). MSESUM 
REAL*8  RHO(26,^6),M^E(5600),DSEED,ALPHA(26,26),BETA(26,26),H(5,5) 
REAL*8  SCALE, YSUM  • 
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c 

C     SET  THE  NUMBER  OF  WHITE  NOISE  REALIZATIONS  USED  TO  CALCULATE 
C  THE  AVERAGE  RHO  MATRIX 

TRIAL  =  25 
C 
C     DEFINE  MODEL  PARAMETERS 


C 


N  =  2 
MN  =  N  * 
MNP1  =  MN  +  1 


C     DEFINE  H(Z)  COEFFICIENTS 
C     THE  AR  PARAMETERS  H(2,l)=0.6  AND  H(l,2)=. 08  CORRESPOND  TO 
C        GAMMA(2)=-. 55  AND  GAMMA(3)=-.08 
DO  7  1=1,  N  v  ' 

DO  6  J=1,N 

H(1,J)  =  0. 

6  CONTINUE 

7  CONTINUE 

C     H(l.l)  =0.5 

H{2,1)  =  0.6 
C     H(tJ)  =  0.2 

Hn,2)  =  0.08 
C     H(2,^)  =  0.2 
C     H(2  3)  =  0.5 
C     H(3,l)  =  0.2 
C 
C     INITIALIZE  MSESUM,  MSEMAX  AND  RHOSUM  MATRIX 

DO  2  I=1,MNP1 
DO  1  J=1AMNP1 

RHOSUM(I,J)  =  0. 

1  CONTINUE 

2  CONTINUE 

MSESUM  =  0. 

MSEMAX  =  0. 

XMAX  =  0. 

YMAX  =  0. 
C     DEFINE  SEED  VALUESA  SATURATION  LIMITX  RANGE  OF  UNIFORM  NOISE., 
C     STD  DEV  OF  GAUSS  NOISE,  AND  NUMBER  OF  POINTS  IN  TIME  SEQUENCES 

IY  =  1354 

DSEED  =  1243073. 5D0 

SAT  =  0.  7 

RANGE  =  1.73205 

STDDEV  =1.0 

NUMPTS  =  5000 

PLTPTS  =  5000 

GRAFN  =  FLOATf NUMPTS  +  1) 

GRAFP  =  FLOAT ( PLTPTS  +  1) 
C 
C    WRITE  HEADER  FOR  UNIFORM  NOISE  INPUT 

WRITE(8.4)  RANGE 

4  FORMATCTf /INPUT  WHITE  UNIFORM  NOISE  HAS  ZERO  MEAN  AND  RANGE  OF 
*+/-  >F6.3) 

WRITEt 8.5V  TRIAL, NUMPTS. I Y 

5  FORMATfTI/RHO  IS  AVERAGED  OVER  ',13,'  TRIALS  OF  ',15,'  POINTS. 
*  INITIAL  SEED=r,I6) 

C 

C    WRITE  HEADER  FOR  GUASSIAN  NOISE  INPUT 
C    WRITE(8+4)  STDDEV 

C4    FORMAT(tl,' INPUT  WHITE  GAUSS  NOISE  HAS  ZERO  MEAN  AND  STD  DEV  OF1 
C    *F6.3) 

C    WRITE(8^5)  TRIALANUMPTS.DSEED 

C5    FORMAT(tl,'RHO  1$  AVERAGED  OVER  ',13,'  TRIALS  OF  ',15,'  POINTS. 
C    *  INITIAL  SEED=t,F10. 1) 
C 

C     PRINT  H  MATRIX 
WRITE(8,8) 

8  FORMATCTf/Y(K)=X(K)-(TENSOR  PRODUCT  H*Y)  WHERE  H(I,J)  =  ') 

DO  10  1=1, N 

WRITE(8,9)  (H(I,J),J=1,N) 
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9  F0RMAT(T10,5(F6.3,4X)) 

10  CONTINUE 

C 

p     **************************************************************** 

C     RUN  TRIALS  FOR  DIFFERENT  SEED  VALUES  TO  GET  AVERAGE  RHO  MATRIX 
p  ****************************** 

DO  50  L=l, TRIAL 

DSEED  =  DSEED/DFLOAT(L) 

IY  =  IY/L 
C     SELECT  EITHER  A  GUASSIAN  OR  UNIFORM  INPUT  NOISE  SEQUENCE 

DO  11  K=1.NUMPTS 
C       THE  INPUT  NOISE  IS  UNIFORM  ON  (-RANGE, RANGE)  WITH 
C       MEAN  =  0  AND  VARIANCE  =((2*RANGE)**2)/12) 

X(K)  =  2.0  *  (URAND(IY7  -  .51  *   RANGE 
C       THE  INPUT  NOISE  IS  GAUSSIAN,  MEAN=0  AND  VARIANCE=STDDEV**2 
C       X(K)  =  GGNQF(DSEED)  *  STDDEV 

11  CONTINUE 
C 

Y(l)  =  DBLE(X(1))  -  H(l.l) 

D02l5=K23LN0Mpfr  "  (  (  '  ]  +   H(2'1}*Y(1)  +  H(3,1)*Y(1)*Y(1)) 

YSUM  =  0. 
DO  14  1=1, N 

DO  13  J=1,N 

IM1  =  1-1 

JM1  =  J-l 

YSUM=H(I,J)*C00RD(Y(K-1),IM1)*C00RD(Y(K-2),JM1)+YSUM 

13  CONTINUE 

14  CONTINUE 

Y(K)  =  DBLE(X(K))  -  YSUM 

15  CONTINUE 
C 

C     ***CALL  SUBROUTINES**** 

C     DETERMINE  AUTOCORRELATION  MATRIX  FOR  Y  SEQUENCE 

CALL  NLCLAT(Y.NUMPTS.R.N) 
C     DETERMINE  REFLECTION  FACTORS  FROM  AUTOCORRELATION  MATRIX 

CALL  SCHURC RHO, R, ALPHA, BETA, MNP1) 
C    ADD  TOGETHER  TH£  R^HO  MATRICES^  FROM  EACH  TRIAL 
DO  45  I=1,MNP1 
DO  40  J=1.MNP1 

RHO$UM(I,J)  =  RHO(I,J)  +  RHOSUM(I,J) 
40      CONTINUE 
45    CONTINUE 
50    CONTINUE 

C     CALCULATE  AVERAGE  RHO  MATRIX  AND  TRUNCATE  TO  TWO  DECIMALS 
DO  54  I=1,MNP1 
DO  53  J=1,MNP1 

RHO(l.J)  =  RHOSUM(I,J)/DFLOAT(TRIAL) 
RHO(I  J)  =  DINT(RHO(I,J)*100. )/100. 

53  CONTINUE 

54  CONTINUE 


C 

55    FORMAT l        LJ)  =') 


WRITE(8V55) 

^f , '  RHO(I .. 

DO  60^1  =  l.MNPl 


WRITEC8,58)CRH0(I,J),J=1,MNP1) 
58      F0RMAT(T1,10F6.2) 
60    CONTINUE 
C 
C 

C     ***  NORMALIZE  THE  LATTICE  INPUT  SIGNALS  AND  PASS  THE  *** 
C         NOISE  GENERATED  DATA  THROUGH  THE  LATTICE  FILTER. 

CALL  NORMS(Y,N,NUMPTS,NORM) 

SCALE  =  1.0/NORM(1) 

CALL  NLLAT( Y , RHO , N , NUMPTS , NORM , XHAT) 

C     ***EXAMINE  THE  WHITENING  EFFECT  OF  THE  LATTICE  FILTER  BY  ***** 
C       CALCULATING  THE  MEAN-SQARE  ERROR  BETWEEN  THE  INPUT  WHITE 
C       NOISE  AND  THE  FORWARD  ERROR  SIGNAL  OUTPUT  FROM  THE  TOP 
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c 


ROW  OF  THE  LATTICE. 
DO  63  K=1,NUMPTS 


MSESUNf  =  fDBLEfXfK))-XHAT(K))*(DBLE(X(K))-XHAT(K))  +  MSESUM 
MSEfK)  =DSQRTfMSESUM/DFLOAT(K)) 
MSEPLTfKl  =  SNGLfMSEfK)) 


IFQISEPLT(K).GT.WSEMAX)  MSEMAX  =  MSEPLT(K) 
63    CONTINUE 
C 

C     NOW  THAT  THE  INVERSE  FILTER  HAS  BEEN  EVALUATED,  DO  THE 
C     DECONVOLUTION  PROBLEM.   FIRST,  DEFINE  THE  "UNKNOWN"  INPUT 
C     SEQUENCE  X(K)  AND  THEN  GENERATE  THE  SYSTEM  OUTPUT  Y(K). 


DO  85  K=1,PLTPTS 

IF(K.  LE.  250)   XfK)  =  FLOAT(K)/250. 

IFfK.  LE.  750.  AND.  K.  GT.  250]  XfK)  =  2.  0  -   FLOATf K)/250. 
IFfK.  LE.  1250.  AND.  K.  GT.  750)  XfK)  =  FLOAT(K1/250.    -4.0 
IFfK.  LE.  1750.  AND.  K.  GT.  1250)   XfK)  =  1.0-   FLOATf K-1250)/250. 
IFfK.  LE.  2000.  AND.  K.  GT.  1750)   XfK)  =  FLOAT(K-1750)/250.    -   1.0 
IFCK.  LE.  2400.  AND.  K.  GT.  2000)   XfK)  =  0.5 


•     IFI 

K.  LE 

IF 

IK.  LE 

85 

CONTINUl 

C 

IFfK.  LE.  2700.  AND.  K.GT.  2400). XfK)  =  -0.5 
IFfK.  LE.  2900.  AND.  K.GT.  2700)   XfK)  =  0.5 
IFfK.  LE.  3000.  AND.  K.GT.  2900)   XfK)  =  -0.5 


3500.  AND.  K.  GT.  3000)XfK)=  0.  7*SIN(0.  0126*FLOAT( K-3001)) 
5000.  AND.  K.GT.  3500)X(K)=  0.  9*SIN(0.  0021*FLOAT(  K-3501) ) 


Yfl)  =  DBLE(Xfl))  -  Hf 1.1) 

Yf2)  =  DBLEfXf2J)  -  (H(i,i)  +  H(2,1)*Y(1)  +  H(3,1)*Y(1)*Y(1)) 
DO  90  K=3,PLTPT$ 
YSUM  =  0. 
DO  89  1=1 .N 

DO  38  J=1,N 

IM1  =  1-1 
JM1  =  J-l 
YSUM=H(I,J)*C00RD(Y(K-1),IM1)*C00RD(Y(K-2),JM1)+YSUM 

88  CONTINUE 

89  CONTINUE 

Y(K)  =  DBLE(X(K))  -  YSUM 

90  CONTINUE 
C 

C     NOW  PASS  THE  SYSTEM  OUTPUT  DATA  THROUGH  THE  LATTICE  FILTER 
C     EMBEDDED  WITH  THE  PREVIOUSLY  CALCULATED  REFLECTION  FACTORS 
C     RHOfI  J).  THE  FORWARD  ERROR  SIGNAL  OUT  OF  THE  TOP  ROW  OF  THE 
C     LATTICE  IS  XHATfKK  AND  SHOULD  APPROXIMATE  THE  DESIRED  SIGNAL 
C     X(J<)  AFTER  IT  IS  RENORMALIZED  AND  SCALED. 

CALL  NORMS(Y,N,PLTPTS,NORM) 

CALL  NLLATfY'RHO,N,PLtPTS,NORM,XHAT) 

DO  98  K=l.PLtPTS  v 

XHATfKf)  =  SCALE  *  f XHAT(K)*NORM(l)) 
XHPLOT(K)  =  SNGLfXHAT(K)) 


YPLOTfK)  =  SNGL(YfK))  /  /  SN 

IF(ABS(XfK)].GT.XMAX)  XMAX  =  ABS(XfK))  t  x 
IFfABS(XHPLOTfK)).GT.XMAX)  XMAX  =  ABSfXHPLOTf K) 
IFf.ABS(YPLOT(K)).GT.  YMAX)  YMAX  =  ABS(YPLOT(K)) 


98  CONTINUE 

XMIN  =  -1.0  *  XMAX 
YMIN  =  -1.0  *  YMAX 
DO  99  K=l,5000 

NINDEtf(K)  =  FLOAT(K-l) 

99  CONTINUE 
C 

C     ****  DISSPLA  PLOTTING  ROUTINES  **** 

C 

C     PLOT  SYSTEM  INPUT  AND  LATTICE  OUTPUT 

C     CALL  TEK618 

CALL  COMPRS 

CALL  RESETf 'ALL') 

CALL  PAGE(8.  0,6.5) 
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I 

i'l'''.' 


; 


CALL  XINTAX 

CALL  AREA2D(6.  75,5.0] 

CALL  XNAMEpTIME  SAMPLES!1  ,100) 

CALL  YNAMEfX(K)   XHAT(K)$'    100) 


CALL  YNAME(/X(K)$\100) 
CALL  CROSS 

CALL  GRAFfO. ,500. ,GRAFP,XMIN, *SCALE' ,XMAX) 
CALL  LINESP   (2.  0) 
CALL  LINES( 'X(Ki$' ,IPAKA1) 
CALL  LINES(.IXHAT(K)$'  ,1PM, 2) 
CALL  LEGLIN 
CALL  DOT 

CALL  CURVE(NINDEX,X,PLTPTS,0) 
CALL  RESETf'DOT1) 

CALL  CURVE(NINDEX,XHPLOTAPLTPTS,0) 
CALL  LEGEND  (IPAK,2,5. 0,6. 5) 
CALL  ENDPL(O) 


C 

C     PLOT  Y 


AREA2D(6.  75,5.0) 
XNAMEPTIME  SAMPI 
.  YNAME('Y(K)$r,10Cv 
CALL  HEADINC'YCK)  =  X(K)  -  0. 2*Y(K-1)**2$' ,100,1.  5.1) 

CALL  HEADIN(7(K)  =  X(K)-(0.  6*Y[ K-1>0.T)8^Y(KL2))V'  100.1.  5,1) 
CALL  HEADIN(hY(K)=X(K)-(0.  I*Y(K-1)+  0. 5*Y(K-1)*Y(K-2J**2 }$' , 


CALL  . 

CALL  XNAMEf'TIME  SAMPLES$' ,100) 

CALL  YNAME('Y(K)$r,100j 


C 

C    *100. 1.5,1) 

CALL  CROSS 

CALL  GRAF(0. ,500. ,GRAFP,YMIN' SCALE' ,YMAX) 

CALL  CURVECNjNDEX  YPLOT  PLTPtS.O) 

CALL  ENDPL(O) 
'  «»            C 
*ti                              C     PLOT  MSE 

CALL  AREA2D(6.  75,5.0) 

CALL  XNAME(VTIME  SAMPLES?' ,100) 

,;::  CALL  YNAMEf'MSEfK)!1 .100] 

CALL  GRAFtD. ,500. ,SrAfN,0. , 'SCALE' ,MSEMAX) 

CALL  CURVEfNlNDEX  MSEPLt,NUMPTS,0) 

CALL  ENDPL(O) 

CALL  DONEPL 
C 

STOP 

END 
p*************  *************************  ********************************* 

C 

C     SUBROUTINE  SCHUR 

C 

C     CALCULATES  THE  REFLECTION  FACTORS  FROM  THE  CORRELATION  MATRIX 

C 

C     WRITTEN  29  APRIL  1985  BY  P.J.  LENK  (REF  12:  P.  174) 

p*********************************************************************** 

SUBROUTINE  SCHUR(RHO. R, ALPHA, BETA. N) 

REAL*8  RHO(26,26),R(26;26),ALPHA(26,26),BETA(26,26),RNORM,T 

C     INITIALIZE  THE  ALPHA  AND  BETA  ARRAYS 
C 

DO  10  I  =  1,N 
DO  5  J  =  l.N 

ALPHA(I\J)  =  R(I,J)/DSQRT(R(I,I)) 

BETA(t.J)  =  ALPHA(1,J) 

RHO(I,J)  '=  0.0 
5       CONTINUE 

C       WRITE(8,7](ALPHA(I,J),J=1,N) 
C7       F0RMAT(5(2X,E12. 5))  ''     * 
10    CONTINUE 
C 
C     BEGIN  CALCULATING  THE  REFLECTION  FACTORS 


C 


DO  50  J  =  2,N 
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NJ1  =  N  -  J  +  1 
DO  40  I  =  1,NJ1 
JI1  =  J  +  I  -  1 
IP1  =  I  +  1 

RH0(I,JI1)  =  ALPHA(I.JI1)/BETA(IP1,JI1) 
RNORM  =  DSQRT(1.0  -  ftHO(I ,JI1)*RH0( I ,JI1)) 
DO  30  K  =  I.N 

T  =  ALPHA(I,K) 

ALPHAfl.K)  =  (ALPHA(I,K)-RH0(I,JI1)*BETA(IP1,K))/RN0RM 
BETA(I,K)  =  (BETA(IP1  K)-RHO(I  JI1)*T)/RN0RM 
30        CONTINUE 
40      CONTINUE 

C        WRITEt8,42)J.((ALPHA(I,K),K<L,N),I=l,N) 
C42      FORMAT<72X'l3:4(2X.E12.5)1 
C       WRITE(8,42)J,((BETA(I,K),K=1,N),I=1,N) 
50    CONTINUE 
C 
C 

RETURN 
END 

p*********************************************************************** 

c 

C  FUNCTION  URAND 

C 

C  TAKEN  FROM  "COMPUTER  METHODS  FOR  MATHEMATICAL  COMPUTATIONS"  BY 

C  G.  E.  FORSYTHE,  M.  A.  MALCOLM,  AND  C.  B.  MOLER 

C 

p********  ********************************************  ******************* 

C 

REAL  FUNCTION  URAND(IY) 

INTEGER  IA.IC,ITWO.M2,M,MIC 

DOUBLE  PRECISION  HALFM 

REAL  S 

DOUBLE  PRECISION  DATAN,DSQRT 

DATA  M2/0/JTWO/2/ 

IF(M2.  NE.  0)  GO  TO  20 
C 

C     IF  FIRST  ENTRY,  COMPUTE  MACHINE  INTEGER  WORD  LENGTH 
C 

M=l 
10  M2=M 

M=ITWO*M2 

IF(M.  GT.  M2)  GO  TO  10 

HALFM=M2 
C 

C     COMPUT  MULTIPLIER  AND  INCREMENT  FOR  LINEAR  CONGRUENTIAL  METHOD 
C 

IA=8*IDINT(HALFM*DATAN( 1.  D0)/8.  D0)+5 

IC=2*IDINT(HALFM*(. 5D0-DSQRT(3.  D0)/6.  D0))+1 

MIC=(M2-IC)+M2 
C 

C     S  IS  THE  SCALE  FACTOR  FOR  CONVERTING  TO  FLOATING  POINT 
C 

S=. 5/HALFM 
C 

C     COMPUTE  NEXT  RANDOM  NUMBER 
C 

20  IY=IY*IA 
C 

C    THE  FOLLOWING  STATEMENT  IS  FOR  COMPUTERS  WHICH  DO  NOT  ALLOW 
C     INTEGER  OVERFLOW  ON  ADDITION 


C 


IF(IY.GT.MIC)  IY=(IY-M2)-M2 
IY=IY+IC 


C 

C    THE  FOLLOWING  IS  FOR  COMPUTERS  FOR  WHICH  THE  WORD  LENGTH 

C     FOR  ADDITION  IS  GREATER  THAN  FOR  MULTIPLICATION 

C 
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IF(IY/2.GT.M2)IY=(IY-M2)-M2 

C    THE  FOLLOWING  STATEMENT  IS  FOR  COMPUTERS  WHERE  INTEGER  OVERFLOW 

C    AFFECTS  THE  SIGN  BIT 

C 

IF(IY. LT.  0}IY=(IY+M2)+M2 

URAND=FLOAT(IY)*S 

RETURN 

END 
p*******  **************************************************************** 

C 

C  SUBROUTINE  NLCLAT                                    * 

C 

C  THIS  SUBROUTINE  PRODUCES  A  CORRELATION  MATRIX  FROM  NONLINEAR   * 

C  TIME  SEQUENCE  IN  AN  ORDER  WHICH  IS  COMPATIBLE  WITH  SUBROUTINE  * 

C  SCHUR.                                              * 

C  * 

C  WRITTEN  7  MAY  1985  BY  P.J.  LENK  (REF  12: P.  181)             * 

I"'  C  * 

p*********************************************************************** 

SUBROUTINE  NLCLAT  (Y.IYS.R.N) 
REAL*8  Y(5000),R(26,26),$UM,VEC(26) 

C    DEFINE  CONSTANTS 
C 

MN  =  N*N 
MNP1  =  MN  +  1 
IYSM2  =  IYS  -  2 
FIYSM2  =  FLOAT(IYSM2) 
!  C 

4  C     INITIALIZE  R  MATRIX  TO  ZERO 


:j  c 


DO  20  I  =  1.MNP1 
DO  10  J  =  1.MNP1 
R(I,J)  =  0.0 
10      CONTINUE 
20    CONTINUE 
C 

C     BEGIN  OUTER  LOOP 
C 

DO  80  I  =  3, IYS 
IR  =  1 

VEC(IR)  =  Y(I) 
DO  50  MP1  =  1,N 
MO  =  MP1  -  1 
LLIM  =  2*MP1  -  1 
DO  40  L  =  l.LLIM 
LO  =  L  -  1 
II  =  MO 
Jl  =  LO/2 

IF  (MOD(L0,2).EQ.  0)  GO  TO  30 
II  =  Jl 
Jl  =  MO 
30  IR  =  IR  +  1 

VEC(IR)  =  C00RD(Y(I-1),I1)*C00RD(Y(I-2),J1) 
40        CONTINUE 
50      CONTINUE 
C 

C     CALCULATE  THE  CORRELATIONS 
C 

DO  70  J  =  1.MNP1 
DO  60  K  =  J.MNP1 

R(J,K)  =  R(J,K)  +  VEC(J)*VEC(K) 
C  WRITE(6'12)VEC(J),VEC(K),R(J,K) 

C12  FORMAT(3(2X,E12. 5)) 

60         CONTINUE 
70     CONTINUE 
80    CONTINUE 
C 
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C     DIVIDE  BY  THE  NUMBER  OF  DATA  ELEMENTS  CONSIDERED 

DO  100  J  =  1.MNP1 

DO  90  K  =  J,MNP1 

R(J,K)  =  R(J,K)/FIYSM2 
90      CONTINUE 
100   CONTINUE 
C 
C     FILL  IN  THE  SYMMETRIC  HALF  OF  CORRELATION  MATRIX 

DO  120  I  =  2,MNP1 
IM1  =1-1 
DO  110  J  =  1.IM1 
R(I,J)  =  R(J,I) 
110      CONTINUE 
120   CONTINUE 
C 
C 

RETURN 
END 
C 
C 

r************************* *********************************  ************* 

C 

C     FUNCTION  COORD 

C 

C     GENERATES  OUTPUT  OF  RANDOM  FUNCTION 

C     CREATED  23  AUG  84  (REF  12: P.  187) 

C 
r*********************************************************************** 

DOUBLE  PRECISION  FUNCTION  COORD(XJ) 
C     USE  SIMPLE  POWER  SERIES  TYPE  POLYNOMIALS 
C 

Y  =  1.0 


IF  (LEO.  0)  GO  TO  30 
Y  =  X**I 


30    COORD  =  DBLE(Y) 

RETURN 

END 
C 
C 

r************************************ *******  ************************** 

C 

C     SUBROUTINE  NLLAT 

C 

C     THIS  SUBROUTINE  IMPLEMENTS  THE  NONLINEAR  LATTICE  FILTER  USING 

C     PREVIOUSLY  CALCULATED  REFLECTION  FACTORS. 

C 

C     WRITTEN  30  APRIL  86 

C 

r*********************************************************************** 


SUBROUTINE  NLLATfY^ RHO, N, NUMPTS, NORM, XHAT) 
:**VARIABLE  DEFINITIONS*^* 


C 
C 

C     INFWD(I)  =  FORWARD  ERROR  INPUT  INTO  THE  I'TH  ROW  OF  THE  LATTICE 
C     INBKDJn.=  BACKWABD. 

C 


OUTFWD(I)=  FORWARD  ERROR  OUTPUT  FROM  THE  I'TH  ROW  OF  THE  LATTICE 
OUTBKDCI)=  BACKWARD  "     "      "   ll    "   "   ''   " 


C  Y  =  INPUT  DATA  VECTOR 

C  RHO  =  REFLECTION  FACTOR  MATRIX 

C  NORM  =  VECTOR  OF  NORMS  OF  THE  LATTICE  INPUT  TERMS 

C  XHAT  =  OUTPUT  DATA  VECTOR;  IT  IS  THE  FORWARD  ERROR  SIGNAL  FROM 

C  THE  LAST  STAGE  OF  THE  FIRST  ROW 

C  NUMPTS  =  NUMBER  OF  POINTS  IN  THE  INPUT/OUTPUT  SEQUENCES 

C  N  =  DIMENSION  OF  SQUARE  Y  DATA  MATRIX 

C  MNP1  =  DIMENSION  OF  THE  RHO,  NORM,  AND  INPUT/OUTPUT  ARRAYS 

C 

C  ***VARIABLE  DECLARATIONS**** 

INTEGER  N,MN, NUMPTS, LAST, MNP1 
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REAL*8  Y(5000).RHO(26,26).NORMf26),XHATf5000),INFWD(26) 
REAL*8  INBKD(2^),0UTFWD(26),0UTBKD(26),RN0RM 

MN  =  N  *  N 
MNP1  =  MN  +  1 
C     INITIALIZE  THE  INPUT  AND  OUTPUT  ARRAYS  TO  ZERO 
DO  5  1=1.26 

INFWOYl)  =  0. 

INBKD(I)  =  0. 

OUTFWD(I)  =  0. 

OUTBKD(I)  =  0. 
5     CONTINUE 
C 
C     DETERMINE  THE  LATTICE  INPUTS  FOR  EACH  TIME  K. 

DO  60  K=1,NUMPTS 
C       CLEAR  INPUTS  FOR  THIS  TIME  ITERATION 
DO  10  I=1,MNP1 

INFWD(I)  =  0. 

INBKD(I)  =  0. 

10  CONTINUE 

IFCK.EQ.  1)  GO  TO  15 
IF(K.  EQ.  2)  GO  TO  20 
C       SET  LATTICE  INPUTS  FOR  CASE  WHEN  K>=3 
IR  =  1 

INFWD(IR)  =  Y(K)/NORM(IR) 
DO  13  MP1  =  1,N 

MO  =  MPi  -  1 
LLIM  =  2*MP1  -  1 
DO  12  L  =  l.LLIM 
LO  =  L  -  1 
II  =  MO 
Jl  =  LO/2 

IF  (MOD(L0,2).EQ.  0)  GO  TO  11 
II  =  Jl 
Jl  =  MO 

11  IR  =  IR  +  1 

INFWD(IR)=C00RD(Y(K-1),I1)*C00RD(Y(K-2),J1)/N0RM(IR) 

12  CONTINUE 

13  CONTINUE 

C  INFWD(l)  =  Y(K)/N0RM(1) 

C  INFWDC2)  =  l./N0RM(2) 

C  INFWD(3)  =  Y(K-1)/N0RM(3' 

C  INFWD(4)  =  Y(K-2)/N0RM(4' 

C  INFWD(5)  =  Y(K-lrY(K-2)/N0RM( 

C  INFWD(6)  =  Y(K-lrY(K-l)/NORM( 

C  INFWD(7)  =  Y(K-2)*Y(K-2l/N0RMl . 

C  INFWD(8)  =  YfK-l)*Y(K-l1*Y(K-2);N0RM(8) 

C  INFWD(9]  =  YCK-l)*Y(.K-2)*Yf  K-2)/N0RM(9) 

C  INFWDC10)  =  Y(K-i)*Y(K-i)*Y(K-2)*Y(K-2)/NORM(10) 
DO  14  J=1,MN 

INBKD(J)  =  INFWD(J+1) 

14  CONTINUE 

GO  TO  25 
C 
C       SET  INPUTS  FOR  K=l  CASE 

15  INFWD(l)  =  Y(1)/N0RM(1) 

"'1(2) 


C 


1. 


INFWD(2)  =  1.7N0RM( 
INBKD(l)  =  INFWD(2) 
GO  TO  25 


C       SET  INPUTS  FOR  K=2  CASE 
INFWD(l)  =  Y(2)/N0RM(1) 
INFWDC2)  =  l./NORMm 

INFWD(3)  =  Y(1)/N0RM(3) 


INFWD(6)  =  Ym*Y(l)/N0RM(6) 
INFWD(ll)  =  Y(1)*Y(1)*Y(1}/N0RMC11) 
INFWDf 18)  =  Yfl1*Y(l)*Y(l)*Y(l)/N0R 


INBKD(l)  =  INFWD(2 
INBKD(2;  =  INFWD(3; 


Y(1)/N0RM(18) 
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INBKD(5)  =  INFWD(6) 
INBKD(IO)  =  INFWD(ll) 
INBKD(17)  =  INFWD(18) 

C     IMPLEMENT  THE  J'TH  UPPER  DIAGONAL  OF  THE  LATTICE.  MOVE  ALONG  THE 
C     DIAGONAL  BY  STARTING  AT  ROW  L=l  AND  MOVING  DOWN  TO  ROW  L=LAST. 
25      DO  50  J=1,MN 

L/ST  =  MNP1  -  J 

IFfj.EO.  1)  GO  TO  40 
C       IF  NOT  ON  THE  FIRST  UPPER  DIAGONAL  (I.E.  J  NOT  =  1)  THEN  UPDATE 
C       THE  INPUTS  FROM  THE  OUTPUTS  OF  THE  PREVIOUS  DIAGONAL 

DO  35  L=1,LAST 

INFWD(L)  =  OUTFWD(L) 

INBKD(L)  =  0UTBKD(L+1) 
35  CONTINUE 

C  DO  THE  LATTICE  CALCULATIONS 

40  DO  45  L=1,LAST 

tfNORM  =  DSQRT(1.0  -  RHO(L,J+L)*RHO(L,J+L)) 
OUTFWD(L)  =  (INFWD(L)  -  RHO(L,J+L)*INBKD( L))/RNORM 
OUTBKD(L)  =  (INBKD(L)  -  RHO(L  J+L)*INFWD(L))/RNORM 
45  CONTINUE 

50      CONTINUE 
C      REMOVE  NORMALIZATION  FACTOR  FROM  THE  OUTPUT  OF  THE  TOP  ROW 

XHAT(K)  =  OUTFWD(l) 
60    CONTINUE 
RETURN 
END 
C 
C 
C 
C 
C 

r* ********************************************************************* 

c 

C     SUBROUTINE  NORMS 

C 

C    THIS  SUBROUTINE  CALCULATES  THE  NORMS  OF  THE  INPUTS 

C    TO  THE  NONLINEAR  LATTICE. 

C 

C    WRITTEN  30  APRIL  86 

C 

r *********************************************************************** 

SUBROUTINE  NORMS(Y,N,NUMPTS,NORM) 

C 

INTEGER  NUMPTS,NAMN,MNP1,I,K,IR,LLIM,L0,M0,MP1,I1,J1,J,L 

REAL*8  VEC(26);n6RM(26),Y(£000) 

MN  =  N  *  N 

MNP1  =  MN  +  1 
C     INITIALIZE  VECTORS  TO  ZERO 

DO  10  K  =  1.MNP1 
VEC(K)  =  0. 
NORM(K)  =  0. 
10    CONTINUE 
C 

C     SINCE  TIME  INDEX  STARTS  AT  1=3^  INITIALIZE  NORMS  OF  INPUTS  WHICH 
C     ARE  POWERS  OF  Y(I-l)  TO  ACCOUNT  FOR  TIMES  1=1  AND  1=2. 

NORM(l)  =  Y(l)  *  Y(l)  +  Y(2)  *  Y(2) 

N0RM(2)  =  1.0+1.0 

N0RM(3)  =  Ym  *  Y(l) 

NORM(6)  =  N0RM(3)  *  NORM(3) 

NORM? 11)  =  NORM(6)  *  NORM(3) 

N0RM(18)  =  NORM(6)  *  NORM(6) 

C     FOR  TIME  U    FORM  A  Y  DATA  VECTOR  WHOSE  COMPONENTS  MATCH  THE 
C     INPUTS  TO  THE  NONLINEAR  LATTICE. 
DO  80  I  =  3,NUMPTS 
IR  =  1 
VEC(IR)  =  Y(I) 
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DO  50  MP1  =  1,N 

MO  =  MPi  -  1 
LLIM  =  2*MP1  -  1 
DO  40  L  =  l.LLIM 
LO  =  L  -  1 
II  =  MO 
Jl  =  LO/2 

IF  (MOD(L0,2).EQ.  0)  GO  TO  30 
II  =  Jl 
Jl  =  MO 
30  IR  =  IR  +  1 

VEC(IR)  =  C00RD(Y(I-1),I1)*C00RD(Y(I-2),J1) 
40  CONTINUE 

50      CONTINUE 
C 
C       CALCULATE  THE  NORMS 

DO  60  K=1,MNP1 
C  IF(K.  Eb.2)   GO  TO  60 

NORM(K)  =  NORM(K)  +  VEC(K)*VEC(K) 
60      CONTINUE 
80    CONTINUE 
C 

WRITE(8,83) 
83    FORMATCTJ'K1  J10,'NORM(K)') 
DO  85  K=l,MNf>l 

NORMCIQ  =  DSQRT(NORM(K)/DFLOAT(NUMPTS)) 
WRITECS,86)  K,NORM(K) 
86      F0RMAT(T1,I3,T8,E12.5) 
85    CONTINUE 
RETURN 
END 
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